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NUMERICAL  MODEL  OF  TWO-DIMENSIONAL  BEACH  CHANGE 


PART  I :  INTRODUCTION 

/ 

1.  The  objective  of  this  report  is  to- present' a  two-dimensional  (or 
schematized  three-dimensional)  numerical  model  of  beach  change  due  to 
breaking  of  short-period  waves.  The  model  development  aimed  at  reproducing 
changes  in  macro-scale  features  of  the  profile,  such  as  bars  and  berms, 
especially  in  the  vicinity  of  coastal  structures.  As  a  first  step  towards 
a  general  numerical  model  for  describing  the  influence  of  structures  on 
beaches,  the  evolution  of  a  beach  around  an  impermeable  groin  or  jetty  was 
simulated. 

2.  A  review  was  made  of  the  literature  on  beach  change  modeling. 
Beach  evolution  models  can  be  classified  into  three  broad  types  as: 

a .  Contour  line  models,  in  which  the  profile  is  schematized  by 
contour  lines  which  monotonically  increase  in  depth  with  distance 
offshore . 

b .  Two-dimensional  (2D)  and  three-dimensional  (3D)  models,  in  which 
microscale  and  mesoscale  features  of  the  physical  processes  (waves, 
currents,  sand  transport)  are  described. 

c .  Cross-shore  transport  (profile  change)  models,  which  simulate 
cross-shore  processes  (beach  profile  change),  but  neglect  longshore 
processes . 

3.  Contour  line  models.  The  shoreline  change  simulation  model 
GENESIS  (Hanson  1987)  defines  the  state  of  the  art  of  1-line  models,  but  it 
does  not  include  cross-shore  sand  transport.  The  zero-depth  contour,  i.e., 
shoreline,  is  taken  as  the  one  line  characterizing  beach  change  in  this 
type  of  model.  One-line  models  have  been  widely  used  in  engineering 
studies,  and  n-line  models  (which  simulate  the  change  in  several  beach 
contours)  have  also  been  developed  to  predict  long-term  (order  of  years) 
shoreline  change  (Perlin  and  Dean  1983).  However,  n-line  models  are 
expensive  to  run  and  have  seen  only  restricted  use  in  engineering  applica¬ 
tions,  probably  because  they  give  only  slightly  more  information  than  thc- 
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1-line  model,  but  at  greatly  increased  computational  cost.  Cross-shore 
sand  transport  formulations  used  in  n-line  models  are  also  greatly  simp¬ 
lified  and  not  well  verified. 

4.  2D  and  3D  models.  There  are  various  levels  of  2D  and  3D  models, 
but  all  solve  the  equation  of  mass  conservation  in  two  horizontal  coor¬ 
dinates.  They  mainly  differ  by  the  sophistication  of  the  wave,  current, 
and  sediment/sand  transport  relationships  used,  and  may  either  compute 
through  the  vertical  or  use  depth- averaged  quantities.  Although  descrip¬ 
tions  of  several  3D  and  quasi- 3D  nearshore  change  models  have  appeared  in 
the  literature,  applications  have  been  limited  to  highly  funded  projects  or 
pure  research  testing,  with  little  verification  against  data  from  real 
beaches.  Because  these  models  are  computationally  intensive,  thev  are 
mainly  applicable  to  projects  of  relatively  small  spatial  extent  i less  than 
approximately  1000  m  in  longshore  extent)  and  short  time  period.,  .  less  that: 
approximately  1  month).  They  also  require  extensive  verification  .d  a 
high  level  of  expertise  when  implemented. 

5.  Cross-shore  (profile  change)  models.  Two  profile  change  modt  1  s 
are  presently  available  for  engineering  use.  One  is  the  model  of  Kriebel 
and  Dean  (1985)  which  treats  the  profile  schematically  and  is  limited  to 
simulation  of  erosional  events.  The  other  model  was  developed  by  Larson 
and  Kraus  (1989).  This  model  describes  the  formation  and  growth  of  bars 
and,  to  a  lesser  extent,  berms,  and  it  can  thus  treat  a  time  sequence  of 
erosion  and  accretion  events, 

6.  On  the  basis  of  the  literature  review  and  the  objective  to 
develop  a  rapid  2D-model  (or  schematized  3D)  applicable  to  projects  of  both 
small  and  large  spatial  extent,  it  was  decided  to  construct  a  2D  grided 
model  which  would  combine  the  simplicity  and  reliability  of  a  shoreline 
change  model  with  a  more  sophisticated  description  of  profile  change.  The 
longshore  transport  rate  calculation  was  expanded  as  compared  to  a  tradi¬ 
tional  shoreline  change  model  to  describe  the  distribution  of  longshore 
sand  transport  across  the  surf  zone  (instead  of  only  a  total  rate  as  used 
in  1-line  models).  The  profile  change  model  is  run  on  shore-normal  lines 
to  compute  cross- shore  sand  transport  rates.  Thus,  in  contrast  to  previou¬ 
sly  developed  2D  and  3D  models,  the  present  model,  named  "3DBEACH,"  for 
schematized  3-Dimensional  BEAch  CHange  model,  is  not  employing  a  complete 
2D  wave  and  current  calculations,  but  locally  decouples  longshore  and 
cross-shore  processes  and  treat  them  individually  as  1-dimensional. 
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7.  The  full  model  3DBEACH  consists  in  principal  of  the  following 
five  modules: 

a.  Wave  propagation  and  decay  module. 

b.  Longshore  current  distribution  module. 

c.  Longshore  transport  module. 

d.  Cross -shore  transport  module. 

e.  Bottom  contour  orientation  module. 

f.  Bottom  topography  change, module . 

8.  In  the  following  chapters  is  first  an  overview  given  of  the 
different  parts  of  the  beach  change  model.  This  involves  a  brief  discus¬ 
sion  of  the  theoretical  foundation  of  each  module  and  the  assumptions 
behind  the  calculations.  Next  chapter  presents  the  numerical  solution 
scheme  and  the  discretization  of  the  governing  equations.  An  short  intro¬ 
duction  is  then  given  to  the  computer  program  as  regards  input  requirements 
and  program  operation.  Finally,  some  sample  calculations  are  shown  to 
illustrate  the  capabiliti-  ,  of  the  model. 
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PART  II:  BEACH  CHANGE  MODEL  OVERVIEW 


9.  The  numerical  model  consists  of  various  modules  which  perform 
calculations  independently  but  are  connected  through  in-  and  output  of  wave 
and  sediment  transport  characteristics.  An  approach  is  used  which  decoupl¬ 
es  individual  profile  lines  when  determining  the  wave  field,  cross-shore 
transport,  and  longshore  transport.  The  interaction  between  profile  lines 
are  represented  through  changes  in  the  bottom  contours  and  conservation  of 
sand.  The  purpose  of  decoupling  the  profiles  in  the  calculations  is  to 
decrease  computational  requirements  since  a  normal  simulation  comprehends 
many  grid  points  and  long  simulation  times.  Solving  the  two-dimensional, 
coupled  wave  equation  in  space,  even  in  its  simplest  form,  put  severe 
restrictions  on  length  and  time  step,  especially  if  steep  slopes  occur 
which  is  the  case  on  the  seaward  side  of  a  bar. 

10.  At  each  time  step  is  the  wave  height  distribution  determined 
along  every  profile  line  including  the  processes  of  shoaling,  refraction, 
and  energy  dissipation  due  to  breaking.  From  the  wave  height  calculations 
are  the  cross-shore  transport  rate  determined  using  transport  relationships 
developed  by  Larson  and  Kraus  (1989).  A  rapid  scheme  for  calculating  the 
longshore  current  (Kraus  and  Larson  1990)  is  used  which  forms  the  basis  for 
determining  the  longshore  transport  together  with  the  wave  height  calcula¬ 
tions  (Kraus,  Gingerich,  and  Dean  1989).  From  the  cross-shore  and  long¬ 
shore  transport  rate  distribution  along  each  profile  line  the  bottom  topog¬ 
raphy  change  is  calculated  from  the  mass  conservation  equation. 

11.  In  the  following  a  brief  description  is  given  of  the  main 
modules  of  the  numerical  model,  that  is,  modules  for  calculation  of  wave 
height,  longshore  current,  cross-shore  transport  rate,  longshore  transport 
rate,  bottom  contour  orientation,  and  bottom  topography  change. 

Wave  Module 


12.  The  wave  height  distribution  is  calculated  along  each  profile 
line  using  linear  wave  theory.  Inside  the  surf  zone,  after  wave  breaking, 
the  breaker  decay  model  by  Dally,  Dean,  and  Dalrymple  (1989)  is  used  to 
calculate  the  wave  height  distribution.  The  point  of  incipient  wave  break¬ 
ing  is  determined  from  an  empirical  criterion  based  on  the  surf  similarity 
parameter  (Battjes  1979)  expressed  in  terms  of  the  deepwater  wave  steepness 
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and  the  local  beach  slope  seaward  of  the  breakpoint  evaluated  over  1/3  of 
the  wavelength. 

13.  An  input  wave  height,  wave  period,  and  incident  wave  angle  must 
be  given  to  the  model  either  at  deepwater  conditions  or  at  a  specific  depth 
corresponding  to  the  gage  depth.  From  this  information  is  the  wave  condi¬ 
tions  calculated  at  the  most  seaward  point  of  the  grid  taking  into  account 
shoaling  and  refraction.  Energy  dissipation  in  the  viscous  layer  along  the 
bottom  is  neglected  to  save  execution  time  since  the  effect  on  the  wave 
height  is  negligible  especially  in  comparison  to  the  amount  of  energy 
dissipation  during  wave  braking  in  the  surf  zone. 

14.  The  calculation  of  the  wave  height  distribution  across -shore 
starts  in  the  most  seaward  calculation  cell  and  proceeds  shoreward  cell  bv 
cell.  The  local  wave  height,  wavelength,  and  wave  angle  with  respect  to 
the  bottom  contours  are  determined  from  Snell's  law  and  the  conservation  of 
energy  flux.  When  wave  breaking  occurs  energy  dissipation  due-  to  increased 
turbulence  in  the  water  column  is  added  in  the  energy  flux  equation.  The 
main  orientation  of  the  bottom  contours  along  a  profile  line  is  determined 
from  a  special  module  described  later  in  this  report. 

15.  The  energy  dissipation  in  the  surf  zone  is  assumed  to  be  propor¬ 
tional  to  the  excess  energy  flux  over  a  stable  energy  flux,  which  is  a 
function  of  the  water  depth  (Horikawa  and  Kuo  1966) .  Wave  properties  in 
the  surf  zone  are  calculated  using  linear  wave  theory  w’ith  no  shallow  water 
approximations.  The  wave  height  calculation  is  carried  out  to  a  water- 
depth  which  is  less  than  the  depth  where  the  surf  zone  ends  in  accordance 
with  the  definition  of  the  different  hydrodynamic  regions  in  the  nearshore. 
The  depth  at  the  shoreward  end  of  the  surf  zone  is  roughly  equal  to  the 
seaward  limit  of  the  backrush . 

16.  The  wave  module  allows  wave  reformation  to  take  place  leading  to 
several  breakpoints  and  thus  surf  zones.  W'ave  reformation  occurs  if  the 
local  wave  height  becomes  lower  than  the  stable  wave  height,  implying  no 
wave  energy  dissipation.  When  the  energy  dissipation  ceases,  shoaling  and 
refraction  will  become  dominant  until  wave  breaking  occur  again  and  dis¬ 
sipation  is  initiated  once  more. 

17.  Generalizing  the  breaker  decay  model  by  Dally,  Dean,  and  Dal- 
rymple  (1984)  the  two-dimensional  equation  for  conservation  of  energy  flux 
incorporating  energy  dissipation  due  to  wave  breaking  may  be  written 
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o  .  co s0)/3x  +  S(F  sinS)/3y  -  K/h  (F  -  Fg)  (1; 

where 

F  -  wave  energy  flux  (Nm/m/sec) 

6  -  wave  angle  with  respect  to  the  bottom  contours 
x  -  cross-shore  coordinate  positive  in  the  seaward  direction  (m) 
y  -  alongshore  coordinate  (m) 
k  =  wave  decay  coefficient 
h  -  water  depth  (m) 

Fg  =  stable  wave  energy  flux  (Nm/m/sec) 

The  wave  energy  flux  is  given  as 

F  -  E  Cg  (2) 

where 

E  =  wave  energy  density  (Nm/m2) 

Cg  -  wave  group  velocity  (m/sec) 

The  wave  energy  density  is  written  using  linear  wave  theory 

E  -  1/8  pgH2  (3) 

where 

H  -  wave  height  (m) 

p  =  density  of  water  (kg/m2) 

g  -  acceleration  of  gravity  (m/sec  ) 

Decoupling  of  the  profile  lines,  equivalent  to  assuming  that  the  wave 
conditions  are  locally  uniform  alongshore  (or  varies  slowly)  and  that  the 
bottom  contours  are  locally  straight  and  parallel  (or  has  a  small 
gradient),  implies  that  Equation  1  may  be  reduced  to 

d(F  cos# )/dx  -  x/h  (F  -  Fs)  (4) 

18,  The  wave  group  velocity  is  related  to  the  velocity  of  an  in¬ 
dividual  wave  C  through  a  factor  n  which  is  a  function  of  the  water  depth 
and  the  wavelength  L  (or  wave  period  T) 

Cg  -  n  C  (5) 
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where 


n  -  1/2  (1  +  (2irh/L)/sinh  (2wh/L))  (6) 

The  speed  of  a  wave  is  determined  through  the  dispersion  relationship 

C  -  CQ  tanh  (2mh/L)  (7) 

where 

C0  =  gT/2jr  (£) 

CQ  is  the  wave  speed  at  deepwater  conditions.  Equation  7  is  an  implicit 
equation  which  has  to  be  solved  numerically  since  C  contains  a  dependence 
on  the  wavelength,  that  is  C  =  L/T. 

19.  The  wave  decay  coefficient  controls  the  rate  of  energy  dissipa¬ 
tion  whereas  the  stable  energy  flux  determines  the  amount  of  energy  dis¬ 
sipation  necessary  in  order  for  stable  conditions  to  occur  once  breaking  is 
initiated.  In  this  respect  stable  wave  conditions  refers  to  a  state  when 
energy  dissipation  due  "o  wave  breaking  ceases  and  thus  waves  reform.  The 
stable  energy  flux  may  be  expressed 

Fs  -  Es  Cg  (9, 

where 

O 

Es  -  stable  wave  energy  density  (Nm/m  ) 

The  stable  wave  energy  flux  corresponds  to  a  stable  wave  height  which  is  a 
function  of  the  water  depth 

Hs  -  Th  (10) 

where 

T  -  stable  wave  height  coefficient 

20.  Two  empirical  coefficients  enter  in  the  breaker  decay  model 
namely  k  and  T.  However,  the  values  of  these  coefficients  seem  to  vary 
negligibly  over  a  wide  range  of  conditions,  both  in  the  laboratory  and  on 
natural  beaches.  Dally,  Dean,  and  Dalrymple  (1984)  recommended  the  values 
of  x  -  0.15  and  T  -  0.40  to  be  used  as  an  average  for  r  varying  bottom 
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slope  based  on  laboratory  and  prototype  -  scale  flume  data.  Ebersole  (1987) 
evaluated  the  performance  of  the  breaker  decay  model  using  field  measure¬ 
ments  and  found  good  agreement  for  engineering  applications.  He  also 
tested  the  breaker  decay  model  developed  by  Svendsen  (1984)  but  found  it  to 
be  inferior  in  terms  of  minimizing  the  root-mean-square  error. 

21.  Larson  and  Kraus  (1989)  used  prototype  -  scale  data  from  Kajima  et 
al .  (1983)  obtained  in  a  large  wave  tank  to  verify  the  applicability  of  the 
breaker  decay  model  for  a  barred  beach  profile.  The  values  obtained  for 
the  two  empirical  parameters  were  much  in  agreement  with  values  obtained  by 
the  aforementioned  authors.  A  significant  dependence  on  beach  slope  for 
the  parameters  was  noted  but  no  specific  relationship  was  developed. 

22.  Kraus  and  Larson  (1989)  evaluated  the  performance  of  the  breaker 
decay  model  for  waves  approaching  the  beach  at  an  angle  with  respect  to  the 
bottom  contours  both  for  laboratory  (Visser  1982)  and  field  data  (Wu, 
Thornton,  and  Guza  1983).  In  these  cases  no  calibration  of  the  model  was 
carried  out  with  respect  to  the  parameters  but  ic  -  0.15  and  F  -  0.40  were 
used.  In  all  cases  did  the  breaker  decay  model  produce  satisfactorily 
agreement.  For  the  field  data  a  statistical  approach  was  used  where  wave 
heights  were  randomly  picked  from  a  Rayleigh  distribution  for  a  large 
number  of  waves  and  the  correspondingly  calculated  wave  height  distribu¬ 
tions  were  averaged  to  obtain  a  representative  distribution. 

23.  As  waves  propagate  onshore  a  flux  of  momentum  (radiation  stress) 
arises  which  causes  displacement  of  the  mean  water  elevation  if  changes  in 
this  flux  occurs  due  to  shoaling  or  breaking.  Outside  the  surf  zone  shoal¬ 
ing  produces  an  increase  in  wave  height  which  causes  a  corresponding  in¬ 
crease  in  the  momentum  flux.  This  flux  increase  is  balanced  by  lowering  of 
the  mean  water  elevation  called  set-down  (Longuet-Higgins  and  Stewart 
1963).  Inside  the  surf  zone,  as  waves  continue  to  break  and  decrease  in 
wave  height,  the  momentum  flux  decreases  and  an  increase  in  mean  water 
elevation  occurs  known  as  set-up. 

24.  The  displacement  of  the  mean  water  surface  (set-up  or  set-down) 
may  be  determined  from  the  momentum  equation 


dSxx/dx  -  -p g(h  +  rj)  drj/dx 


(11) 


Sxx  -  radiation  stress  in  the  direction  of  the  waves  (N/m) 
i)  -  wave  set-up  or  set-down  (m) 

The  radiation  stress  is  given  by  linear  wave  theory  for  an  arbitrary  wave 
angle  of  incidence 

Sxx  -  1/8  pgH2  (n  (cos2  6  +  1)  -  1/2)  (12) 

25.  The  energy  flux  conservation  equation  (Equation  4),  the  wave 
dispersion  relationship  (Equation  7),  and  the  momentum  equation  (Equation 
11)  is  not  sufficient  to  specify  the  wave  height  distribution  across-shore 
if  the  waves  have  an  angle  of  incidence  with  respect  to  the  bottom  con¬ 
tours.  In  this  case  an  additional  equation  has  to  be  solved  which  expres¬ 
ses  the  conservation  of  wave  number.  The  wave  number  conservation  equation 
is  written  in  the  general  form,  not  taking  into  account  diffraction 

5(sin«/L )/3x  -  3(cos#/L)/3y  -  0  (131 

26.  For  straight  and  parallel  bottom  contours,  implying  no  (or  mild- 
alongshore  variation  in  wave  and  bottom  characteristics,  Equation  13  may  be 
simplified  to  yield 

d(sin0/L)/dx  -  0  (14) 

If  Equation  14  is  integrated  Snell's  law  is  obtained  stating  that  the 
quantity  sinS/L  is  a  constant  across-shore.  Equations  4,  7,  11.  and  14 
allows  the  wave  height  distribution  to  be  determined  once  the  boundary 
conditions  have  been  established. 

27.  As  mentioned  previously  the  wave  characteristics,  such  as  wave 
height,  period,  and  incident  angle,  must  be  specified  at  the  seaward  end  of 
each  profile.  The  end  of  the  grid  should  be  located  seaward  of  the  break 
point  at  all  times.  From  the  wave  properties  it  is  possible  to  determine 
analytically  the  set-down  at  the  most  seaward  point  of  the  profiles  (Lon- 
guet-Higgins  and  Stewart  1962) 

rj  -  irH2/(4L  sinh(4?rh/L) )  (13) 

28.  Also,  a  condition  for  the  initiation  of  breaking  is  needed 
indicating  when  energy  dissipation  due  to  breaking  starts  occurring  and  the 
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wave  decay  coefficient  becomes  different  from  zero.  Breaking  occurs  when 
the  ratio  between  the  wave  height  and  the  water  depth  exceeds  a  certain 
value  which  is  a  function  of  the  surf  similarity  parameter  expressed  in 
terms  of  the  deepwater  wave  steepness  and  profile  slope  prior  to  breaking 

7  -  1.14  f0'21  (16) 

where 

7  -  breaker  ratio  (-H^/h^,  where  index  b  denotes  breaking  conditions) 

6  -  surf  similarity  parameter 
The  surf  similarity  parameter  is  given  by 

?  -  tan^/jH \JT0  (18) 

where 

fi  -  local  beach  slope  seaward  of  break  point 

29.  Equation  16  was  derived  by  Larson  and  Kraus  (1^89)  from  data 
given  by  Kajima  et  al.  (1983)  but  very  similar  relationships  were  previous¬ 
ly  reported  in  Singamsetti  and  Wind  (1980)  and  Sunamura  (1980b).  Larson 
and  Kraus  (1989)  used  data  from  prototype  -  scale  experiments  whereas  the  two 
latter  studies  were  based  on  small-scale  laboratory  experiments.  Thus,  in 
Equation  16  scale  effects  seem  to  be  absent  although  it  should  be  observed 
that  only  monochromatic  waves  were  used  in  the  tests  with  probably  neglig¬ 
ible  wave-by-wave  interaction. 

30.  The  experimental  data  by  Kajima  et  al.  (1983)  encompass  cases 
with  a  maximum  wave  period  of  12  sec  although  most  of  the  cases  have  per¬ 
iods  in  the  range  3-9  sec.  The  limited  data  material  for  longer  period 
waves,  often  encountered  in  the  field,  makes  it  necessary  to  use  Equation 
16  with  care.  Also,  since  the  relationship  was  determined  mainly  from 
barred  profiles  exposed  to  monochromatic  waves,  the  profile  slope  seaward 
of  the  break  point  was  in  general  steeper  than  what  is  found  in  the  field. 
Thus,  Equation  16  tended  to  underestimate  the  breaker  ratio  in  the  field 
for  gentle  slopes  (Larson  and  Kraus  1989). 

31.  The  energy  dissipation  due  to  wave  breaking  is  needed  in  the 
prediction  of  the  cross-shore  sand  transport  rate.  As  the  wave  height 
decays  after  breaking  turbulent  motion  is  induced  in  the  surf  zone.  The 
turbulent  kinetic  energy  is  converted  into  heat  through  transfer  of  energy 
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from  larger  to  smaller  eddies.  In  the  macro-scale  approach  of  the  beach 
profile  model  it  is  satisfactory  to  derive  energy  dissipation  from  wave 
height  decay  without  further  resolving  the  internal  structure  of  this 
process.  However,  the  difference  between  the  production  of  turbulent 
kinetic  energy  and  the  corresponding  dissipation  with  reference  to  the 
time-scale  (Roelvink  and  Stive  1989)  is  observed  as  well  as  the  process  of 
energy  reordering  (Svendsen,  Madsen,  and  Buhr  Hansen  1978). 

32.  The  wave  energy  dissipation  D  per  unit  water  volume  may  be 
expressed  as 

D  -  1/h  d(F  costf)/dx  (19) 

Using  Equation  4  the  following  relationship  for  D  is  obtained 

D  -  k/ h2  (F  -  Fs)  (20) 

Once  the  wave  height  distribution  is  calculated,  the  energy  dissipation  per 
unit  water  volume  is  determined  explicitly  from  Equation  20  in  the  surf 
zone.  Outside  the  surf  zone  k  is  regarded  to  be  zero  and  no  energy  dis¬ 
sipation  takes  place  since  bottom  friction  is  neglected. 

Longshore  Current  Module 

33.  When  the  wave  height  distribution  has  been  determined  from  the 
wave  module  the  longshore  current  distribution  along  each  profile  line  is 
calculated.  A  simplified  approach  involving  a  rapid  numerical  solution 
scheme  presented  by  Kraus  and  Larson  (1990)  is  used.  The  governing  equa¬ 
tion  is  the  momentum  equation  in  the  longshore  direction  (the  momentum 
equation  in  the  cross-shore  direction  is  already  used  in  the  wave  module  to 
compute  the  set-up)  which  is  a  balance  between  the  change  in  the  alongshore 
radiation  stress  directed  onshore,  the  bottom  friction,  and  the  lateral 
mixing  due  to  turbulent  diffusion. 

34.  The  governing  equation  may  be  written 

A  d^v/dx^  +  dA/dx  dv/dx  -  Bv  ■=  1/p  dSv„/dx  (21) 

where 

v  =  longshore  current  velocity  (m/sec) 
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Sv„  -  alongshore  radiation  stress  directed  onshore  (N/m) 

The  coefficients  A  and  B  in  Equation  21  is  (Kraus  and  Larson  1990) 

A  -  €  (h  +  v)  (22) 

B  -  -2/ir  CfUjj,  (1  +  sin2  6)  (23) 

where 

2 

€  -  lateral  mixing  coefficient  (m  /sec) 

C£  -  bottom  friction  coefficient 

Uj^  -  maximum  bottom  orbital  velocity  (m/sec) 

35.  The  lateral  mixing  coefficient  is  proportional  to  the  local 
maximum  orbital  velocity  and  the  local  wave  height  according  to 

c  -  Av^H  (29) 

where  the  proportionality  coefficient  is  denoted  as  A  and  has  a  typical 
value  of  1.0.  The  bottom  friction  coefficient  normally  has  a  value  of 
about  0.01  for  a  sandy  beach. 

36.  The  alongshore  radiation  stress  directed  onshore  decreases  in 
the  surf  zone  as  the  wave  height  decreases  due  to  breaking  and  provides  the 
main  driving  force  for  the  longshore  current.  This  stress  is  expressed  as 

Sxy  -  1/16  pgH2  n  sin  26  (25) 

The  lateral  mixing  term  in  Equation  21  flattens  the  longshore  current 
distribution  and  produces  a  velocity  outside  the  break  point  where  no 
change  in  radiation  stress  occurs. 

37.  The  longshore  current  module  were  extensively  tested  both  again¬ 
st  laboratory  and  field  data  to  estimate  coefficient  values,  and  to  verify 
the  predictive  capabilities  of  the  module.  The  same  data  sets  as  were  used 
for  the  verification  of  the  wave  module  were  applied  (Visser  1982,  Wu, 
Thornton,  and  Guza  1985)  together  with  data  from  Kraus  and  Sasaki  (1979) 
and  Ebersole  and  Dalrymple  (1980).  Ebersoie  and  Dalrymple  (1980)  presented 
a  hypothetical  calculation  example  for  a  barred  beach  which  provided  a  test 
for  the  current  module  regarding  numerical  stability  for  a  complex  bottom 
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topography  and  also  a  possibility  to  compare  with  a  more  sophisticated 
model  of  nearshore  circulation. 

38.  The  rapid  and  simplified  longshore  module  used  in  the  beach 
change  model  proved  to  reproduce  both  the  laboratory  and  field  data  in  a 
satisfactory  manner.  Figure  1  displays  a  comparison  between  the  calculated 
longshore  current  distribution  and  the  current  distribution  measured  by 
Visser  (1982)  in  one  of  his  experiments  (case  2).  Also,  the  calculated  and 
measured  wave  height  distribution  is  shown.  The  parameter  values  in  the 
breaker  decay  model  recommended  by  Dally,  Dean,  and  Dalrymple  (1984)  were 
used  without  any  calibration. 


Longshore  Current  Speed  (m/sec)  Wave  Height  (m) 


Figure  1.  Comparison  between  calculated  and  measured  wave  height  respectively 
longshore  current  distribution  using  data  from  Visser  (1982) 

39.  For  the  field  data,  two  different  techniques  were  employed  to 
determine  the  longshore  current  distribution.  In  the  first  case  a  repres¬ 
entative  wave  height,  the  root-mean-square  wave  height,  was  used  in  the 
calculations.  In  the  other  case  a  large  number  of  wave  heights  were 
randomly  picked  from  a  Rayleigh  distribution  and  the  current  distribution 
corresponding  to  each  wave  height  was  determined.  To  obtain  a  representa¬ 
tive  current  distribution,  the  average  of  all  calculated  distributions  were 
taken.  In  comparison  with  the  field  data  a  somewhat  better  agreement  was 
obtained  with  the  first  method  if  the  breaker  ratio  was  put  to  0.4. 
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However,  this  value  is  to  be  regarded  as  a  statistical  description  of  the 
conditions  at  an  inferred  break  point,  which  encompasses  both  broken  and 
unbroken  waves. 

40.  Accordingly,  the  second  technique  was  preferred  in  the  simula¬ 
tions  since  it  clearly  distinguish  between  broken  and  unbroken  waves  which 
is  of  fundamental  importance  when  determining  the  sediment  transport  rate. 
The  simulations  showed  that  random  forcing  conditions  in  the  field  could  be 
reproduced  by  superimposing  the  effect  of  a  number  of  monochromatic  waves. 
For  the  case  with  a  complex  bottom  topography  (Ebersole  and  Dalrymple  1980) 
the  model  proved  to  be  inherently  stable  and  show  satisfactory  qualitative 
agreement  with  the  complete  circulation  model. 

Cross-Shore  Sediment  Transport  Module 

41.  The  cross -shore  movement  of  sediment  in  a  surf  zone  is  governed 
by  the  properties  of  the  velocity  field  and  the  sediment  concentration. 
Thus,  a  description  of  the  velocity  and  concentration  field  is  necessary 
not  only  across  the  surf  zone  but  also  through  the  water  column  at  every 
point.  At  our  present  stage  of  knowledge,  such  a  detailed  prediction  is 
not  possible  in  engineering  numerical  models  of  beach  change  to  be  used  on 
natural  beaches.  The  sediment  concentration  in  the  surf  zone  is  closely 
related  to  the  generation  of  turbulent  motion,  which  depends  on  the  wave 
breaking.  Accordingly,  it  is  plausible  that  the  sediment  concentration, 
and  thus  the  amount  of  material  available  for  transport,  is  closely  related 
to  the  wave  energy  dissipation  due  to  wave  breaking. 

42.  In  the  cross-shore  sand  transport  module,  no  prediction  of  the 
cross-shore  velocity  field  is  made  but  a  simple  engineering  approach  is 
used  to  determine  the  transport  rate  distribution.  The  direction  of  the 
transport  is  predicted  from  an  empirical  criterion  derived  from  the  large 
wave  tank  data  and  verified  with  field  data,  whereas  the  magnitude  of 
transport  is  given  by  the  wave  energy  dissipation  per  unit  water  volume. 
This  simplifying  approach  implies  a  transport  rate  distribution  directed 
either  on-  or  offshore  along  the  entire  profile  at  a  specific  instant  in 
time . 

43.  Criterion  for  predicting  the  transport  direction  or  the  overall 
profile  evolution  has  been  suggested  by  a  number  of  authors  (for  a  brief 
review,  see  Larson  and  Kraus  1989).  This  type  of  criterion  normally 
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comprehends  one  parameter  characterizing  the  incident  wave  conditions  and 
one  parameter  involving  some  property  of  the  sediments  (grain  size,  D,  or 
fall  speed,  w) .  In  the  present  work,  the  criterion  developed  by  Larson  and 
Kraus  (1989)  was  used  involving  the  two  parameters  H0/LQ  and  H0/wT  (su¬ 
bscript  o  denotes  deepwater  conditions) ,  which  produced  a  good  distinction 
between  profiles  exhibiting  on-  and  offshore  sand  movement,  both  in  the 
laboratory  and  in  the  field. 

44.  In  Figure  2  is  the  criterion  by  Larson  and  Kraus  (1989)  dis¬ 
played  (the  solid  line)  together  with  data  from  large  wave  tank  experi¬ 
ments.  A  satisfactory  division  is  achieved  described  by  the  equation 

Ho/Lo  ‘  M  (Hq/wT)3  (26) 

where 

H  -  wave  height  (m) 

L  —  wavelength  (m) 
w  -  sand  fall  speed  (m/sec) 

T  —  wave  period  (sec) 

in  which  M  -  0.0070  and  the  subscript  o  refers  to  deepwater  conditions.  If 
the  left  side  of  Equation  26  is  less  (greater)  than  the  right  side,  the 
profile  is  predicted  to  erode  (accrete) . 

45.  In  the  wave  tank  studies  the  choice  of  wave  and  sediment 
properties  to  be  used  in  Equation  26  is  trivial  due  to  constant  wave 
forcing  conditions  and  well-defined  sediment.  However,  in  order  to  apply 
the  equation  to  a  natural  beach  the  suitable  wave  and  sediment  characteris¬ 
tics  have  to  be  found.  Field  data  was  therefore  collected  from  various 
investigations  around  the  world  to  establish  the  proper  statistical 
descriptors  to  be  used  in  Equation  26.  The  data  used  are  summarized  in 
Sunamura  (1980a)  and  Seymour  (1986),  where  the  movement  of  the  deeper 
contours  were  used  to  evaluate  profile  change  and  transport  direction  in 
the  latter  investigation. 

46.  Several  different  statistical  wave  characteristics  (mean  wave 
height  H,  root-mean-square  wave  height  Hrms>  significant  wave  height  Hso) 
were  used  in  order  to  obtain  the  best  division  between  on-  and  offshore 
(accretion/erosion)  sand  movement.  The  mean  wave  height  proved  to  give  the 
best  delineation  using  the  same  coefficient  values  as  for  the  large  wave 
tank  data,  that  is,  the  wave  height  in  Equation  26  should  be  taken  as  the 
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Figure  2.  Criterion  for  predicting  cross-shore  sand  transport  direction  based 
on  monochromatic  laboratory  data 

mean  for  field  application.  Figure  3  displays  in  analogy  to  Figure  2  the 
criterion  (solid  line)  together  with  the  observed  data  from  natural 
beaches.  In  the  case  significant  wave  height  was  given,  the  mean  wave 
height  was  calculated  as  R  «  0.626  Hso  from  assuming  a  Rayleigh  distribu¬ 
tion  to  be  valid  for  the  wave  height.  The  mean  or  peak  period  of  the 
spectra  was  chosen  as  the  representative  period  to  be  used  in  Equation  26 
together  with  the  median  grain  size  of  the  beach  sediment. 

47.  Based  on  the  division  of  the  profile  applied  in  nearshore  wave 
dynamics  and  the  physical  characteristics  of  sediment  transport  under 
various  flow  conditions,  four  different  zones  of  cross-shore  transport 
along  a  beach  profile  were  introduced  by  Larson  and  Kraus  (1989).  These 
zones  are 


Zone  I : 

Zone  II: 

Zone  III: 


From  the  seaward  depth  of  effective  sand  transport  to  the 
break  point  (pre -breaking  zone) 

From  the  break  point  to  the  plunge  point  (breaker 
transition  zone) 

From  the  plunge  point  to  the  point  of  wave  reformation  or 
to  the  swash  zone  (broken  wave  zone) 
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i 
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Figure  3.  Criterion  for  predicting  cross-shore  sand  transport  direction 
applied  to  field  conditions 

Zone  IV:  From  the  shoreward  boundary  of  the  surf  zone  to  the 

shoreward  limit  of  runup  (swash  zone) 

Figure  4  gives  a  definition  sketch  of  the  different  transport  zones. 

48.  The  zone  between  the  break  point  and  the  plunge  point  is 
analogous  to  the  outer  region  (vortex  region  according  to  Hiller  1976) 
proposed  by  Svendsen,  Madsen,  and  Buhr  Hansen  (1978).  A  certain  distance 
is  required  after  breaking  before  the  turbulent  conditions  are  approxima¬ 
tely  uniform  throughout  the  water  column.  From  a  cross-shore  transport 
point  of  view,  assuming  the  transport  rate  being  proportional  to  the  energy 
dissipation  per  unit  water  volume,  the  recognition  of  this  zone  is  impor¬ 
tant  and  physically  motivated.  A  more  detailed  discussion  of  the  different 
transport  zones  are  given  in  Larson  and  Kraus  (1989).  In  this  report  only 

■*’  the  empirically-based  transport  relationships  for  the  zones  derived  from 

the  large  wave  tank  data  will  be  presented.  These  formulas  are  used  in  the 
numerical  model  to  determine  the  transport  rate  distribution  across-shore . 

49.  The  transport  rate  relationships  derived  based  on  physical  con¬ 
siderations  and  analysis  of  the  large  wave  tank  data  may  be  summarized  for 
the  four  different  zones 
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€  -  slope-related  transport  rate  coefficient  (m^/sec) 

h  -  water  depth  (m) 

The  subscripts  b,  p,  s,  and  r  stands  (in  order)  for  quantities  evaluated  at 
the  break  point,  plunge  point,  end  of  the  surf  zone,  and  runup  limit.  Two 
different  spatial  decay  coefficients  are  used  in  Zone  I  and  II  (with  the 
subscripts  1  and  2  respectively)  to  describe  the  decrease  in  transport  rate 
with  distance. 

50.  Empirical  expression  for  the  spatial  decay  coefficients  were 
derived  from  large  wave  tank  data  (Larson  and  Kraus  1989)  .  The  decay 
coefficient  in  Zone  I  was  related  to  the  median  grain  size  and  the  breaking 
wave  height 

*1  =  10.4  (D50/Hb)0-47  (51) 

In  Zone  II  limited  data  suggested  that 

A2  -  0.2  Aj_  (32) 

When  calculating  the  transport  rate  in  Zone  I  and  II  the  transport  rate  is 
first  determined  at  the  plunge  point  from  Equation  29  and  then  the  exponen¬ 
tial  decays  are  applied  seaward  in  respective  zone. 

51.  The  cross-shore  transport  rate  in  zones  of  fully  broken  waves  is 
related  to  the  wave  energy  dissipation  per  unit  water  volume.  This  type  of 
transport  formula  has  previously  been  used  by  Moore  (1982),  Kriebel  (1982, 
1986),  and  Kriebel  and  Dean  (1985).  Analysis  of  large  wave  tank  data 
substantiated  this  type  of  transport  relationship  in  zones  of  broken  waves 
(Larson  1988,  Larson  and  Kraus  1989).  The  transport  relationships  in  the 
other  zones  are  empirical  and  based  on  a  large  amount  of  data  from  wave 
tank  experiments.  For  example,  in  the  swash  zone  the  transport  rate  simply 
decreases  linearly  from  the  end  of  the  surf  zone  to  the  runup  limit.  This 
property  of  the  transport  rate  distribution  appeared  for  many  large  wave 
tank  cases,  both  with  on-  and  offshore  sand  movement,  and  similar  charac¬ 
teristics  have  been  noted  on  natural  beaches  where  the  beach  face  receded 
uniformly  during  erosional  and  accretionary  events  (Seymour  1987). 
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52.  Based  on  the  wave  height  and  longshore  current  calculations  the 
longshore  sand  transport  is  determined  using  the  transport  formula  sug¬ 
gested  by  Kraus,  Gingerich,  and  Dean  (1989).  This  formula  was  developed 
from  measurements  carried  out  during  two  major  field  experiments  at  Duck, 
North  Carolina.  The  transport  rate  is  proportional  to  the  local  wave 
height  and  longshore  current  velocity 

qy  -  k  (VH  -  Cc)  (33) 

where  k  and  Cc  are  empirically  determined  coefficients.  In  Figure  5  is  a 
comparison  shown  from  Kraus,  Gingerich,  and  Dean  (1989)  of  transport  rates 
measured  in  the  field  and  Equation  33,  where  the  wave  height  was  taken  as 
the  root-mean-square  height. 


Figure  5.  Comparison  between  measured  and  predicted  longshore  sand  transport 
rates  (after  Kraus,  Gingerich,  and  Dean  1989) 

53.  Kraus,  Gingerich,  and  Dean  (1989)  also  proposed  more  complex 
equations  for  the  longshore  transport  rate  incorporating  the  spatial 
derivative  of  the  wave  height  and  the  standard  deviation  of  the  current 
velocity.  The  latter  quantity  is  not  possible  to  estimate  but  has  to  be 


measured  in  the  field,  and  thus  can  not  be  used  in  a  predictive  model.  The 
term  containing  dH/dx  was  included  initially  but  omitted  later  since  it 
became  dominant  around  the  break  point  where  the  energy  dissipation  is 
large . 

54.  The  longshore  transport  module  was  extensively  tested  for 
various  hypothetical  wave  and  water  level  conditions  together  with  the  wave 
and  longshore  current  module.  A  maximum  in  transport  rate  normally  occurs 
somewhat  shoreward  of  the  breakpoint  where  the  product  of  the  local  wave 
height  and  longshore  current  has  its  maximum.  Thus,  if  monochromatic  waves 
and  a  constant  water  level  are  employed,  and  beach  change  updrift  an 
impermeable,  infinitely  long  groin  (or  jetty)  is  studied,  the  impoundment 
will  be  largest  close  to  the  break  point.  As  time  elapse  a  bar-like 
formation  will  build  up  at  the  groin  which  eventually  prevents  the  waves 
from  propagating  towards  the  shoreline. 

55.  If  cross -shore  sand  transport  is  taken  into  account  this  build¬ 
up  could  be  reduced  as  material  is  transported  offshore  depending  on  the 
relation  between  the  magnitude  of  longshore  and  cross -shore  transport.  In 
general,  however,  strong  impoundment  at  the  groin  around  the  break  point  is 
to  be  expected  due  to  the  distribution  of  the  longshore  sedimenc  transport. 
This  build-up  normally  does  not  occur  in  the  field  but  is  to  be  attributed 
to  constant  forcing  conditions,  and  the  problem  is  circumvented  by  intro¬ 
ducing  some  variability  in  the  wave  or  water  level  conditions. 

56.  The  longshore  sand  transport  is  assumed  to  occur  along  the 
average  contour  of  a  specific  profile  line.  This  implies  that  the  long¬ 
shore  transport  has  a  component  directed  in  the  cross -shore  direction  when 
reference  is  made  to  the  global  rectangular  coordinate  system.  Thus,  in 
the  beach  change  model  is  a  calculated  longshore  transport  decomposed  with 
respect  to  the  average  bottom  contour  giving  rise  to  a  contribution  to  the 
cross-shore  transport. 

Bottom  Contour  Orientation  Module 

57.  The  bottom  contour  orientation  is  important  in  the  calculation 
of  the  wave  height  distribution  along  each  profile  line.  The  incident 
angle  between  the  wave  crests  and  the  bottom  contours  determines  wave 
refraction  and  also  enters  in  the  wave  energy  conservation  equation. 
Furthermore,  since  the  contour  orientation  is  dependent  upon  the  local 
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depth  conditions  at  several  neighboring  profile  lines,  the  contour  scheme 
introduces  a  coupling  between  the  beach  evolution  along  the  calculation 
grid. 

58.  The  contouring  scheme  in  the  numerical  model  uses  a  linear  plane 
to  approximate  the  depth  surface  around  the  grid  point  where  the  local 
slope  of  the  contour  is  to  be  estimated.  From  the  linear  plane  equation  is 
the  contour  slope  readily  calculated  and  the  average  slope  given  by  two 
planes  is  used  to  obtain  a  slope  estimate  at  the  correct  grid  point 
location  for  the  wave  height  calculations.  The  general  equation  for  a 
linear  plane,  used  to  locally  approximate  the  bottom  topography,  is 

h(x,y)  -  A  +  Bx  +  Cy  (34) 

where  A,  B,  and  C  are  coefficients  given  by  the  depths  in  three  neighboring 
grid  points.  The  local  slope  is  then  given  by 

dy/dx  -  -B/C  (35) 

The  ratio  between  the  coefficients  in  Equation  35  may  be  expressed  in  terms 
of  grid  point  depths  and  distances  as 

B/C  -  (Ay  (h2  -  h1))/(Ax  (2h3  -  hx  -  h2))  (36) 

where 

Ax  -  cross -shore  length  step 
Ay  -  longshore  length  step 

A  definition  sketch  explaining  the  location  of  the  depths  used  in  Equation 
36  is  given  in  Figure  6. 

59.  The  slope  estimate  given  by  Equation  36  is  taken  somewhat  to  the 
right  of  the  actual  grid  point  location  used  in  the  wave  height  calculation 
(see  Part  III).  To  avoid  this  another  slope -estimate  is  derived  from  a 
linear  plane  placed  through  a  different  set  of  grid  point  depths  yielding  a 
left-adjusted  slope  approximation  (see  Figure  6).  The  average  between  the 
two  slope-estimates  provides  a  value  at  the  desired  location. 

60.  The  average  contour  orientation  is  computed  based  on  the  local 
slope  at  each  grid  point  for  every  profile  line  and  used  in  the  wave  height 
calculations.  Initially,  the  local  slope  at  the  grid  p>  int  was  used  but 
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Figure  6.  Definition  sketch  for  determining  local  bottom  contour  slope 

this  approach  caused  difficulties  at  some  grid  points  for  profile  lines 
updrift  a  groin  where  strong  accretion  occurred.  Instead  the  average 
contour  orientation  was  employed  which  proved  to  introduce  more  stability 
in  the  calculations.  The  use  of  an  average  slope  implies  that  the  waves 
respond  to  the  general  characteristics  of  the  bottom  topography  not  to 
local  changes  in  beach  bathymetry,  which  is  in  agreement  with  the  overall 
modeling  approach. 


Beach  Change  Module 


61.  The  change  in  bottom  topography  is  calculated  from  the  mass 
conservation  equation,  which  is  written 


db/dt  -  dqx/dx  +  5qy/5y  (3 

Thus,  once  the  cross-shore  and  longshore  transport  rate  distribution  have 
been  determined,  the  corresponding  beach  change  is  given  by  Fn-'ation  37. 
Through  the  sand  conservation  equation  is  also  the  intera'  i  tween 

neighboring  profile  lines  introduced. 

62.  To  avoid  excess  steepening  of  profile  slopes  during  conditions 
when  bar  growth  is  significant,  under  near  steady  wave  and  water  level 
forcing,  avalanching  is  allowed  to  occur.  The  concept  of  avalanching,  as 
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discussed  by  Allen  (1970),  was  Incorporated  in  the  numerical  model  and  a 
routine  was  included  to  account  for  the  transport  induced  by  slope  failure. 
Allen  (1970)  recognized  two  different  limiting  profile  slopes,  namely  angle 
of  initial  yield  (angle  of  repose)  and  residual  angle  after  shearing.  If 
the  local  slope  exceeds  the  angle  of  initial  yield  (^)  ,  material  is 
redistributed  along  the  slope  and  a  new  stable  slope  is  attained,  known  as 
the  residual  angle  after  shearing  (y6r)  • 

63.  Analysis  carried  out  by  Larson  and  Kraus  (1989)  showed  values  of 
the  maximum  slope,  inferred  as  the  angle  of  initial  yield,  having  a  mean 
value  of  28  deg.  The  average  value  for  the  residual  angle  after  shearing 
was  determined  to  be  22  deg,  however,  the  analysis  was  difficult  to  carry 
out  since  some  profile  steepening  probably  occurred  between  the  avalanching 
and  the  nearest  profile  survey  in  time.  Allen  (1970)  determined  from 
experiments  the  difference  between  ipt  and  0r  (dilatation  angle)  to  be 
in  the  range  10-15  deg  for  sand.  Therefore,  in  the  beach  change  model  is 

i> ^  set  to  28  deg  and  4>r  set  to  18  deg,  that  is,  a  dilatation  angle  of  10 
deg  is  used. 

64.  If  the  local  slope  exceeds  the  material  will  be  redistribu¬ 

ted  in  the  neighboring  cells  so  that  the  slope  everywhere  is  less  or  equal 
V>r.  Computationally  this  is  carried  out  by  checking  the  local  slope  ip  at 
each  time  step  along  the  entire  grid  and  if  xf>  >  Vj  avalanching  is 
initiated.  Two  conditions  determine  the  redistribution  of  material,  (1) 
mass  should  be  conserved,  and  (2)  the  residual  angle  after  shearing  is  \tr. 
In  Figure  7  a  definition  sketch  is  given  where  the  local  slope  between  cell 
1  and  2  exceeds  Depths  before  avalanching  is  denoted  with  h,  whereas 

the  new  depths  after  sand  has  been  redistributed  is  indicated  with  a  prime. 

65.  The  number  of  cells  N  over  which  the  avalanching  is  occurring 
is  not  given  a  priori  but  has  to  be  found  iteratively.  Initially,  it  is 
assumed  that  the  avalanching  takes  place  between  cell  1  and  2  only. 

Material  is  moved  from  cell  1  to  2  in  order  to  obtain  the  slope 

between  the  cells  and  the  new  slope  between  cell  2  and  3  is  checked  if  it 
exceeds  ip{-  If  this  is  the  case,  the  calculations  has  to  be  redone,  now 
incorporating  cell  1-3,  using  the  mass  conservation  equation  and  assigning 
the  slope  Vr  to  all  slopes  between  cells.  This  procedure  is  carried  out 
until  the  slope  between  cell  N  and  N+l  is  less  than  0r. 

66.  For  the  general  case  were  sand  is  to  redistributed  in  N 
calculation  cells  included  in  the  avalanching,  the  mass  conservation 
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Figure  7.  Definition  sketch  for  describing  avalanching  along  the  profile 
equation  is  written 

N 

l  Ah*  -  0  (38) 

i-1 

where 

Ahj.  -  hi  -  hi  (39) 

Observe  that  the  cell  numbering,  denoted  by  i,  starts  from  the  cell  where 
avalanching  is  initiated  and  is  taken  in  the  direction  of  sand 
redistribution.  Also,  the  depth  changes  have  to  be  taken  with  sign.  The 
other  criterion  necessary  to  yield  the  solution  is  that  after  avalanching 
the  slope  between  calculation  cells  has  to  be  ^r.  If  the  length  of  the 
calculation  cell  is  Ax,  then  the  residual  angle  after  shearing  maybe  ex¬ 
pressed  as 

(6r  -  (ht  -  hi+1)/Ax  (40) 

29 


\ 


Since  is  a  constant,  the  difference  in  depth  between  neighboring  cells 

Ah  has  to  be  a  constant  given  by 

Ah  -  dx  (41) 

67.  The  requirement  regarding  the  difference  in  depth  between  cells 
as  given  by  Equation  41  may  be  written  for  N-l  pair  of  cells  included  in 
the  avalanching 

Ah  -  (h-^  +  Ah^)  -  (h2  +  Al^) 

Ah  *  (I12  +  Ah2 )  -  (hj  +  Ahj) 

Ah  -  (ht  +  Ahi)  -  (hi+1  +  Ahi+1)  (42) 

Ah  -  (hfj.]^  +  Ahjj_^)  -  (hfj  +  Ah^) 

This  system  of  equations  may  be  rewritten  by  consecutively  eliminating  the 
first  part  of  the  right-hand  side  of  each  equation  by  adding  the  previous 
equation.  The  result  is 


Ah  -  (h-^  +  Ahj)  -  (h2  +  Ah2> 

2Ah  -  (h-^  +  Ah-^)  -  (hj  +  Ahj) 

iAh  -  (hjL  +  Ahx)  -  (hi+1  +  Ahi+1)  (43) 

( N - 1 ) Ah  -  (h-^  +  Ahj^)  -  (hN  +  AhN) 

68.  From  Equation  43  the  unknown  depth  changes  Ah^  may  be  solved 
for  in  every  equation 


Ah2  _ 

(hx  + 

AhL)  - 

h2  -  Ah 

Ah3  - 

(hi  + 

Ahx)  - 

b3  -  2Ah 

Ahi+i 

-  (hx 

+  Ah;) 

-  hi+i  -  iAh 

AhN  - 

(ht  + 

AhL)  - 

hN  -  (N-l)Ah 

(44) 
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From  Equation  44  it  is  possible  to  compute  all  the  depth  changes  once  the 
depth  change  in  the  first  cell  is  known.  To  determine  the  depth  change  in 
the  first  cell  the  mass  conservation  equation,  Equation  38,  is  used.  The 
sum  of  all  the  changes  in  depth,  taken  with  sign,  has  to  be  zero 

N  N 

£  Ahj  -  Ahx  +  Y  Ahx  —  0  (45) 

i-1  i-2 

Replacing  the  summation  over  the  depth  changes  from  cell  2  to  N  using 
Equation  44  yields 

N  N  N 

Ahx  t  Y  (hl  +  Ah].)  '  I  hi  •  £  Ah  “  0  (46) 

i-2  i-2  i-2 

69.  From  Equation  46  it  is  possible  to  explicitly  solve  for  the 
depth  change  in  cell  1 

N 

Ahx  -  -(N-l)/N  hx  +  (Y  hx)/N  +  (N-l)  Ah/2  (47) 

i-2 

The  depth  changes  in  cell  2  to  N  is  then  given  by 

Ahx  ”  ^1  +  A*^l  ‘  ^i  '  (i'l)  A^  (48) 
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PART  III:  NUMERICAL  SOLUTION  OF  GOVERNING  EQUATIONS 


70.  In  this  chapter  a  brief  description  is  given  of  the  numerical 
techniques  for  solving  the  governing  equations  in  the  beach  change  model. 
The  description  is  primarily  given  for  the  wave  module,  longshore  current 
module,  and  beach  change  module.  Also,  the  boundary  conditions  in  the 
model  are  discussed  and  typical  time  and  length  step  used  in  the  calcula¬ 
tions  are  presented. 

71.  In  Figure  8  is  a  definition  sketch  given  of  the  calculation  grid 
and  where  the  different  quantities  are  taken.  The  wave  height  and  long¬ 
shore  current  computations  are  carried  out  in  between  profile  lines  using 
interpolated  depths  from  the  original  grid  (Original  depths  are  denoted  as 
h  in  Figure  8  whereas  interpolated  depths  between  profile  lines  are  named 
havm  and  have  depending  on  location  in  the  grid.  Note  that  in  the  rest  of 
the  report,  especially  in  the  equations,  depth  is  denoted  as  h  even  if  it 
is  an  interpolated  depth  for  the  sake  of  generality;  the  context  should 
explain  if  an  interpolated  depth  is  used  when  calculations  are  done) . 
Furthermore,  before  employing  the  mass  conservation  equation  the  cross - 
shore  transport  rates  has  to  be  interpolated  in  order  to  obtain  values  at 
the  proper  grid  locations. 

72.  The  index  i  is  used  for  cross-shore  grid  point  numbering  and 
index  j  for  longshore  grid  point  numbering,  The  x-axis  is  positive 
pointing  offshore  and  the  y-axis  is  directed  alongshore  according  to  a 
right-hand  coordinate  system  with  depth  h  positive  below  mean  sea  level. 

Wave  Module 

73.  An  explicit  finite -difference  scheme  is  used  to  solve  the 
equations  for  determining  the  wave  height  distribution  across-shore . 

In  the  model  a  quasi-stationary  approach  is  applied  indicating  that  a  time- 
independent  solution  of  the  wave  height  distribution  is  determined  at  every 
time  step.  The  propagation  of  individual  waves  across-shore  is  not 
described  but  it  is  tacitly  assumed  that  the  model  time  step  is  significan¬ 
tly  larger  than  the  wave  period.  For  irregular  wave  conditions  the  wave 
height  used  at  each  time  step  may  be  regarded  as  a  statistical  measure 
which  represents  the  average  properties  of  the  breaking  waves  during  the 
time  step. 
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Figure  8.  Definition  sketch  of  numerical  grid  for  beach  change  model 

74.  The  numerical  calculation  starts  at  the  end  of  the  grid  and 
proceeds  onshore  for  every  profile  line  through  an  explicit  scheme  where 
quantities  known  at  a  specific  grid  point  are  used  to  determine  correspond¬ 
ing  quantities  at  the  next  grid  point  closer  to  shore.  A  staggered  grid  is 
used  with  different  quantities  taken  at  points  in  the  middle  of  calculation 
cells  (along  a  line  in  between  the  depths  of  the  original  grid)  or  at  the 
boundaries  between  cells  (see  Figure  8).  The  main  quantity  in  the  middle 
of  a  cell  is  the  water  depth  and  these  grid  points  will  be  referred  to  as 
x-points  (see  Figure  8).  At  the  boundaries  of  calculation  cells  the  main 
quantity  is  the  cross-shore  transport  rate  and  the  grid  points  located  here 
are  known  as  #-points  (see  Figure  8). 

75.  From  knowing  all  wave  properties  at  the  first  grid  point  (#- 
point),  the  wave  set-down  is  determined  from  Equation  11.  The  profile 
depth  (havm)  is  given  in  the  middle  of  each  calculation  cell  whereas  the 
set-down  (set-up)  is  calculated  at  #-points.  The  depths  at  the  boundaries 
of  the  cells  (have)  are  given  by  linear  interpolation.  Also,  the  energy 
flux  and  the  radiation  stress  in  the  direction  of  the  waves  is  determined 
from  knowledge  of  the  depth  and  wave  quantities  at  the  first  #-point. 

76.  When  going  from  a  specific  #-point  to  the  next,  wave  refraction 
is  first  determined  if  the  incident  wave  approaches  with  an  angle  to  the 
bottom  contours.  From  Equation  14  the  angle  between  the  wave  crests  and 
the  bottom  contours  at  the  next  #-point  onshore  is  given  by 
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0 i  -  arcsin  (L^/Lj^  sin 


(49) 


Note  that  the  grid  point  numbering  increases  in  the  seaward  direction  since 
the  x-axis  points  offshore,  but  the  calculation  proceeds  from  offshore  in 
the  shoreward  direction.  The  water  depth  at  grid  point  i  in  Equation  49 
is  corrected  with  the  set-down  (set-up)  using  the  value  at  grid  point  i+1 
in  order  to  avoid  an  implicit  scheme.  The  wavelengths  are  calculated  using 
the  dispersion  relationship  (Equation  7)  numerically  solved  by  a  Pade' 
approximation  (Chen  and  Thompson  1985) . 

77.  In  Equation  49  is  9  the  angle  between  the  incident  wave  crests 
and  the  bottom  contours,  whereas  the  bottom  contour  orientation  with 
reference  to  the  coordinate  system  is  denoted  9C  and  the  corresponding  wave 
crest  orientation  denoted  0V.  Thus,  9  -  0V  -  9  c  with  both  angles  being 
positive  counterclockwise  from  the  x-axis.  As  is  pointed  out  in  Figure  8, 
0C  is  taken  in  a  x-point  whereas  9V  is  calculated  in  ^-points. 

78.  Next  step  in  the  calculation  is  to  determine  the  energy  flux  and 
thus  the  wave  height.  Equation  4  is  written  in  difference  form 

Fi  -  (F1+1  (cosfli+1  -  0.5Aci)  +  Ac  F^/cos^  (50) 

where 

Aci  -  k  Ax/(hi  +  iji+1)  (51) 

Ax  -  cross-shore  length  step  (m) 

The  stable  wave  energy  flux  is  determined  from 

Fsi  -  1/8  pg  (r  (hi  +  f?i+1)  )2  (C8i  +  Cgi+1)/2  (52) 

An  average  value  for  the  wave  group  velocities  is  used  since  this  quantity 
is  taken  at  #-points  and  the  stable  wave  energy  flux  is  evaluated  at  x- 
points.  Outside  the  break  point,  k  should  be  set  to  zero,  implying  that 
Aci  is  also  zero,  since  no  energy  dissipation  due  to  breaking  is  occurring. 

79.  When  the  energy  flux  has  been  calculated  at  a  specific  point  the 
corresponding  wave  height  is  determined  from  Equations  2  and  3 
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(53) 


Using  the  wave  height,  the  radiation  stress  is  calculated  from  Equation  12 
and  the  set-down  (set-up)  is  given  from  Equation  11  expressed  in  difference 
form 

'll  "  ''i+l  +  <Sxxi+l  -  sxxi>/^8<hi  +  ''i+l>> 

As  before  rj  is  taken  at  a  #-point  half  a  cell  seaward  of  the  x-point  when 
determining  the  depth  to  avoid  an  implicit  scheme. 

80.  At  every  calculation  step  a  check  is  made  if  wave  breaking 
occurs  according  to  the  criterion  (Equation  16) .  Once  breaking  is  in¬ 
itiated  the  wave  decay  coefficient  is  given  the  value  x  =  0.15  and  energy 
dissipation  is  taking  place.  Also,  if  breaking  has  occurred  a  check  is 
made  if  wave  reformation  is  possible,  that  is,  if  the  wave  energy  flux  is 
lower  than  the  stable  energy  flux.  The  stable  wave  height  coefficient  is 
set  to  T  -  0.40  and  if  there  is  wave  reformation  x  is  again  given  the 
value  zero.  An  arbitrary  number  of  surf  zones  with  intermediate  zones  of 
wave  reformation  may  occur. 

81.  The  wave  height  calculations  are  not  carried  out  for  the  actual 
profiles  but  a  somewhat  smoothed  profile  is  used.  The  motivation  for  this 
is  that  the  response  of  the  wave  to  the  bottom  topography  takes  place  over 
a  certain  distance  suggesting  that  shoaling  and  refraction  does  not  occur 
instantaneously.  For  instance,  on  a  natural  beach  a  small  hump  in  the 
bottom  would  in  general  not  cause  a  wave  to  shoal  up,  which  the  equations 
would  predict  if  no  profile  smoothing  was  applied.  The  number  of  calcula¬ 
tion  cells  over  which  the  smoothing  is  performed  (a  moving  average  scheme) 
is  determined  based  on  the  breaking  wave  height  as  3Hb,  where  the  index  b 
denotes  breaking  conditions.  A  simple  predictive  formula  (Larson  and  Kraus 
1989)  is  used  to  approximately  estimate  the  breaking  wave  height  at  each 
time  step  prior  to  determining  the  wave  height  distribution 

Hb  -  0.525  Hq  (H0/L0)-°-25  (55) 

82.  The  breaker  ratio  is  determined  from  Equation  16  for  the  most 
seaward  break  point.  If  multiple  breaking  occurs  a  constant  breaker  ratio 
7  -  1.0  is  applied  at  the  consecutive  more  shoreward  break  points.  The 
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shoreward  boundary  of  the  calculation  grid  is  a  predetermined  water  depth 
approximately  corresponding  to  the  seaward  limit  of  the  backrush  or  the 
location  of  a  seawall. 


Longshore  Current  Module 

83.  The  longshore  current  distribution  is  determined  along  the  same 
line  as  where  the  wave  height  calculations  are  carried  out.  The  current 
velocity  is  taken  at  x-points  whereas  the  alongshore  radiation  stress 
directed  onshore  is  given  at  #-points  from  the  wave  height  computations.  A 
double -sweep  solution  scheme  is  employed  for  solving  the  tridiagonal  system 
of  equations  that  arises  when  writing  Equation  21  in  difference  form  for 
all  grid  points  along  a  profile  line.  For  grid  point  i  the  following 
equation  is  obtained  written  in  canonical  form 


-AHiVi-1  +  BHivi  -  CHivi+l  " 

(56) 

where 

AH^  -  -Aj/Ax2 

(57) 

BHi  -  -(B1  +  (Ai+1  +  AjVAx2) 

(58) 

CH^  -  -Ai+i/Ax2 

(59) 

^i  "  (sxyi+l  ‘  Sxyi^/Ax 

(60) 

The  quantities  A  and  B  are  defined  in  Equations  22  and  23,  both  taken  at  *- 
points . 

84.  By  introducing  the  two  double-sweep  coefficients  and  F^,  a 
recursive  equation  in  the  current  velocity  may  be  derived 

vi  "  Eivl+1  +  Fi  <61) 

By  substituting  Equation  61,  written  for  index  i  and  i-1,  into  Equation  60, 
the  following  expressions  are  obtained  for  the  double-sweep  coefficients 

Ef  -  CHj/DNi  (62) 

Fi  -  +  RH^/DNj  (63) 
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where 


DNt  -  BHj^  -  AH|Ei_1  (64) 

The  solution  procedure  is  thus  to  determine  all  E^  and  starting  from  i=2 

through  the  whole  grid  with  knowledge  about  the  shoreward  boundary  condi¬ 
tion  at  grid  point  i-1.  The  current  velocities  are  then  calculated  from  a 
backward  sweep  by  using  Equation  61  and  the  seaward  boundary  condition. 

85.  The  seaward  boundary  condition  is  no  current  velocity  at  the  end 
of  the  grid,  whereas  the  shoreward  boundary  condition  is  a  current  velocity 
corresponding  to  the  explicit  solution  of  Equation  21  neglecting  mixing. 

The  current  velocity  at  the  most  shoreward  grid  point  may  be  written 

vj.  -  -1/Bi  dsxy/dx  (65) 

This  is  a  more  realistic  description  of  the  current  conditions  at  the 
shoreline  as  compared  to  assuming  a  zero  velocity  at  the  shoreward  boun¬ 
dary.  From  Equations  61  and  65  the  values  of  the  double-sweep  coefficients 
at  the  shoreward  boundary  is  given  as  E^-0  and  F^-v-^.  Note  that  the  most 
shoreward  grid  point  in  the  longshore  current  calculation,  given  the  number 
1  locally  in  the  current  module,  corresponds  to  a  higher  grid  point  number 
in  the  global  calculation  grid.  The  longshore  current  calculation  proceeds 
shoreward  as  far  as  the  wave  height  computation. 

Beach  Change  Module 

86.  From  the  calculations  carried  out  in  the  wave  and  longshore 
current  module  it  is  possible  to  determine  the  cross-shore  and  longshore 
transport  rate  distribution  along  each  profile  line.  The  cross-shore 
transport  rate  distributions  are  given  by  Equations  27-30  and  the  longshore 
transport  rate  distribution  by  Equation  33.  In  this  section  a  brief 
outline  is  given  of  the  numerical  methods  for  computing  the  respective 
distributions  and  the  discretization  of  the  mass  conservation  equation. 

87.  Cross-shore  transport  distribution.  From  knowledge  of  the 
location  of  the  various  transport  zones  it  is  possible  to  calculate  the 
cross-shore  transport  rate  distribution  from  Equation  27-30.  Several  break 
points  may  occur  along  a  profile  if  wave  reformation  takes  place  leading  to 
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several  zones  of  type  II  and  III.  In  order  to  determine  the  transport  rate 
distribution ,  the  transport  rate  is  first  calculated  in  zones  of  fully 
broken  waves  (Zone  III)  according  to  Equation  29  written  in  difference  form 

qxi  -  K(  (Dj  +  D^p/2  -  Deq  +  e/K  (h£  -  h^p/Ax)  (66) 

The  energy  dissipation  per  unit  water  volume  is  taken  at  x-points  (in  the 
middle  of  a  calculation  cell)  making  it  necessary  to  average  D  when 
calculating  the  transport  rate. 

88.  The  calculation  of  the  transport  rate  in  zones  of  fully  broken 
waves  provides  the  transport  rate  at  the  boundaries  of  Zone  III  from  which 
the  transport  rate  may  be  calculated  in  the  other  zones.  When  the  trans¬ 
port  rate  at  the  plunge  point  and  at  the  end  of  the  surf  zone  is  known, 
Equations  27,  28,  and  30  are  applied  to  completely  determine  the  transport 
rate  distribution. 

89.  Longshore  transport  rate  distribution.  From  knowledge  about  the 
wave  height  and  longshore  current  distribution  the  longshore  transport  rate 
is  calculated  using  Equation  33  written  in  difference  form 

qyi  -  k(Vi(Hi+1  +  Ht)/2  -  Cc)  (67) 

The  calculated  transport  rate  is  regarded  as  being  directed  along  the 
average  contour  orientation  of  a  profile  line  and  decomposed  according  to 
the  global  coordinate  system  before  the  mass  coti__rvation  is  employed. 

90.  Mass  conservation  equation.  The  mass  conservation  equation  is 
solved  explicitly  to  obtain  the  bottom  change  at  every  grid  point.  The 
longshore  transport  rate  is  given  at  the  same  location  as  the  longshore 
current  (see  Figure  8)  may  thus  be  used  directly  in  the  mass  conservation 
equation.  Regarding  the  cross -shore  transport  rate  and  the  contribution 
from  the  longshore  transport  to  the  cross-shore  transport,  interpolation 
has  to  be  carried  out  to  obtain  the  transport  rate  at  the  proper  location. 
The  mass  conservation  in  difference  form  is  expressed  as 

<hij  -  hi,j)/^t  “  (<ixi+i, j  -  qxi.j)/**  +  (qyi.j+i  '  qyi,j)/Ay  (68) 

At  -  time  step  (sec) 

Ay  -  longshore  length  step  (m) 
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where  the  index  k  denotes  a  specific  time  level  during  which  all  transport 
rates  are  taken. 

91.  Since  the  cross-shore  transport  rate  distribution  is  determined 
from  different  relationships  depending  on  the  zone  of  transport,  the 
derivative  will  in  general  not  be  continuous  at  the  boundaries  between  the 
zones.  To  obtain  a  smoother,  more  realistic  transport  distribution  a  three 
point  filter  is  applied  to  the  calculated  cross -shore  transport  rates 
according  to 

9x i  -  9xi-l/^  +  9xi/2  +  <5xi+l/4  <69> 

where  the  prime  denotes  the  smoothed  transport  rate.  The  same  type  of 
expression  is  also  applied  in  the  longshore  direction  to  achieve  a  more 
gentle  variation  in  qy,  especially  in  the  vicinity  of  a  groin  or  jetty. 

92.  The  boundary  conditions  in  the  beach  change  calculation  differs 
somewhat  in  the  cross-shore  and  longshore  direction.  For  each  profile  line 
is  the  cross-shore  boundary  conditions  no  sand  transport  passed  certain 
grid  points  defined  by  the  length  of  the  active  profile.  The  grid  point 
corresponding  to  the  location  of  runup  constitutes  the  shoreward  boundary, 
whereas  the  seaward  boundary  is  determined  by  how  rapid  the  decay  in 
transport  rate  seaward  of  the  break  point  is.  When  the  cross -shore 
transport  rate  has  decreased  to  a  small  predetermined  value  the  calculation 
is  stopped  and  the  current  grid  point  is  regarded  as  the  seaward  boundary. 

93.  Two  different  boundary  conditions  are  possible  to  specify  in  the 
longshore  direction,  to  be  discussed  more  in  Part  IV.  The  transport  rate 
may  be  assumed  to  be  zero  at  a  boundary,  corresponding  to  an  impermeable 
groin,  or  the  boundary  may  be  regarded  as  an  open  coast  with  unaffected 
transport,  that  is,  no  beach  change  at  the  boundary  due  to  longshore 
transport  (pinned  beach).  A  typical  time  step  in  the  model  is  5-20  min,  a 
cross-shore  length  step  of  1-5  m,  and  a  longshore  length  step  of  20-50  m. 
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PART  IV:  ORGANIZATION  OF  COMPUTER  PROGRAM 


94.  This  chapter  describes  the  organization  of  the  computer  code  and 
how  to  operate  the  program.  Details  are  given  regarding  the  input  data 
requirements,  the  external  file  handling,  and  the  output  produced  when 
running  the  program.  The  description  of  the  computer  code  is  aiming  at 
being  extensive  enough  to  serve  as  a  users  manual  for  applying  the  program 
to  beach  change  simulation.  In  Appendix  A  is  the  complete  computer  program 
presented  with  a  large  number  of  comment  statements  to  facilitate  the 
understanding  of  the  various  parts  of  the  code.  Also,  in  Appendix  B  is  an 
example  of  an  input  file  to  the  program  given. 

Computer  Program  Overview 

95.  The  code  is  separated  into  a  number  of  different  subroutines 
performing  specific  tasks  in  the  beach  change  calculations.  In  Figure  9  is 
an  overview  given  of  the  subroutines  and  the  interaction  between  the 
different  routines  during  the  calculations.  The  main  program  routine  is 
called  3DBEACH  and  administrates  the  in-  and  output  of  data.  Initially 
general  information  necessary  for  the  calculation  is  read  in,  to  be 
discussed  in  detail  in  the  next  section,  together  with  the  depth  matrix. 
Also,  3DBEACH  organizes  the  entire  calculation  and  calls  at  every  time  step 
the  subroutines  required  for  computing  the  beach  change. 

96.  The  depth  matrix  is  read  in  by  the  subroutine  DEPIN  with  ten 
depth  values  given  at  every  line  in  the  input  file,  with  respective  profile 
line  appearing  as  columns  in  the  matrix.  WAVIN  obtains  the  wave  and  water 
level  information  at  every  time  step  by  reading  wave  height,  incident  wave 
angle,  wave  period,  and  water  level  with  reference  to  mean  sea  level  (water 
level  positive  above  mean  sea  level).  The  input  wave  conditions  refer  to 
deepwater  conditions  or  a  specific  depth. 

97.  The  wave  height  distribution  for  every  profile  line  is  deter¬ 
mined  in  WAVEHD,  which  uses  the  service  routine  LINWAV  to  calculate  wave 
properties  at  a  specific  depth.  CONPLA  is  also  called  by  WAVEHD  in  order 
to  determine  the  bottom  contour  orientation.  After  the  wave  height 
computations  have  been  carried  the  longshore  current  distribution  is 
determined  through  the  routine  VELOC.  The  cross-shore  and  longshore 
transport  rate  distribution  is  calculated  in  TRANX  respectively  TRANY. 
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Figure  9.  Flow  chart  describing  the  subroutines  in  the  beach  change  computer 
program 


TRANX  uses  the  subroutines  DOMAIN  to  determine  the  locations  of  the 
different  cross-shore  transport  regions,  TRS1GN  to  decide  the  direction  of 
transport,  and  XTRANS  to  calculate  the  transport  rates  at  every  grid  point. 
In  BATHYM  is  the  change  in  bottom  topography  obtained  and  by  calling  the 
routine  MAXANG  is  a  check  performed  throughout  the  grid  if  the  angle  of 
repose  is  exceeded. 


Input  Data  Requirements 


98.  Besides  the  depth  matrix  and  wave  and  water  level  information  at 
every  time  step,  as  previously  discussed,  some  general  computational  data 
have  to  be  supplied  to  run  the  program.  This  data  should  be  entered  into  a 
start  file  (file  handling  is  discussed  in  the  next  section)  which  is  read 
by  the  program  before  any  calculations  are  performed.  In  the  start  file  a 
comment  line  describing  the  input  separates  every  line  of  input  data. 
Appendix  B  shows  the  typical  structure  of  an  input  file  containing  the 
necessary  computational  information. 

99.  At  the  top  of  the  input  file  are  four  lines  of  heading  which  is 
skipped  by  the  program  when  reading  input  data.  The  first  line  of  input 
allows  a  description  of  the  current  simulation  to  be  printed  out  encom¬ 
passing  up  to  70  characters.  Next  line  involves  the  number  of  time  steps, 
numoer  of  grid  points  in  the  cross-shore  direction,  and  number  of  grid 
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points  in  the  longshore  direction.  Thereafter  should  the  time  step  and 
length  step  in  respective  coordinate  direction  be  given.  The  coefficients 
required  for  the  longshore  current  computations  must  be  supplied,  that  is, 
the  bottom  friction  coefficient  and  the  mixing  parameter. 

100.  The  median  grain  size  of  the  beach  is  entered  at  the  next  line 
in  units  of  mm.  The  output  interval,  expressed  in  terms  of  calculation 
loops,  determines  when  output  information,  such  as  bottom  topography  and 
transport  rates,  is  written  to  different  files.  A  depth  has  to  be  supplied 
where  the  offshore  wave  data  were  collected  or  if  deepwater  conditions 
prevail  a  negative  number  should  be  given  on  this  line.  On  the  next  line 
the  coefficients  appearing  in  the  longshore  respective  cross-shore  sand 
transport  relationship  is  entered.  The  final  line  in  the  start  file  refers 
to  the  longshore  boundary  conditions  at  the  first  and  last  profile  line  of 
the  grid.  The  number  0  corresponds  to  an  open  coast  whereas  the  number  1 
implies  an  impermeable  infinitely  long  groin  or  jetty,  blocking  all  sand 
transport . 

101.  For  some  coefficients  input  values  are  hardwired  in  the 
computer  program  or  given  in  data  statements.  In  general,  these  are 
coefficients  with  values  that  are  not  expected  to  change  from  one  simula¬ 
tion  to  another,  but  may  be  assumed  to  be  fairly  constant.  Examples  of 
such  coefficients  are  the  two  parameters  in  the  breaker  decay  model  and  the 
exponential  decay  parameters  in  the  cross-shore  transport  rate  distribution 
(Larson  and  Kraus  1989)  . 


File  Handling  in  the  Program 

102.  A  number  of  sequential  files  with  specific  names  which  handles 
the  in-  and  output  of  data  are  needed  when  operating  the  program.  The 
different  files  are  assigned  logical  names  through  open  statements.  Three 
files  are  needed  for  the  input  of  data  and  three  files  are  used  for  the 
output  from  the  model.  The  file  names  given  in  this  report  refer  to  a  VAX 
computer  system  but  the  names  could  easily  be  modified  for  implementation 
on  any  other  system. 

103.  The  general  input  governing  a  specific  simulation  should  be 
stored  on  a  file  called  START.DAT.  As  discussed  previously,  the  input  data 
is  entered  in  free  format  with  commentary  lines  in  between  lines  with  data 
(see  Appendix  B) .  The  depth  matrix  is  contained  in  a  file  named  DEPTH.DAT, 
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where  each  profile  line  corresponds  to  a  column  in  the  matrix  and  ten 
values  are  given  per  line.  Wave  and  water  level  conditions  are  punched  in 
free  format  in  the  file  WAVES.DAT.  Wave  height,  incident  wave  angle,  wave 
period,  and  water  level  is  given  at  every  time  step  entered  in  free  format 
at  one  line  per  time  step. 

104.  The  output  from  the  program  consists  of  bottom  topography  or 
transport  rates  at  a  chosen  time  increment  throughout  the  simulation 
together  with  the  conditions  at  the  last  time  step.  Other  information  from 
the  program,  such  as  wave  height  and  longshore  current  distribution,  has 
been  prepared  for  printing  to  various  logical  units  but  is  tentatively 
commented  out  to  reduce  the  amount  of  output  from  a  simulation.  The  bottom 
topography  is  written  to  the  file  DEPUT.DAT  in  the  same  format  as  the  input 
depth  matrix.  The  cross- shore  and  longshore  transport  rate  distribution 
along  each  profile  line  are  stored  in  the  files  TRANX.DAT  and  TRANY.DAT 
respectively.  This  is  done  in  the  same  format  as  for  the  depth  files. 
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PART  V:  SAMPLE  CALCULATIONS  WITH  THE  BEACH  CHANGE  MODEL 


105.  In  this  chapter  a  number  a  sample  calculations  with  the  beach 
change  model  are  presented  in  order  to  illustrate  some  of  the  characteris¬ 
tics  of  the  model.  In  all  cases,  the  beach  evolution  updrift  a  groin  was 
simulated  permitting  no  sand  transport  past  the  groin.  At  the  other  end  of 
the  grid  an  open  coast  was  assumed  (pinned  beach)  with  enough  sand  supply 
to  prevent  any  beach  change  due  to  longshore  transport.  Accordingly, 
cross-shore  sand  transport  may  still  occur  at  the  profile  line  located  on 
the  boundaries. 

106.  A  spatially  distorted  grid  was  used  with  a  length  increment  of 
2  m  in  the  cross-shore  direction  and  20  m  in  the  longshore  direction.  The 
time  step  was  set  to  20  min  and  the  total  simulation  time  encompassed  160 
time  steps.  Only  a  limited  grid  was  used  to  reduce  execution  times  on  the 
computer  namely,  100  grid  points  across-shore  and  10  profile  lines  along¬ 
shore.  The  initial  bottom  topography  was  assumed  to  be  uniform  alongshore 
with  a  linear  dune  face  and  a  Dean-type  equilibrium  profile  (Dean  1977) 
from  the  shoreline  and  seaward. 

107.  Two  different  wave  and  water  level  conditions  were  studied 
where  the  wave  height  was  varied  sinusoidally  in  one  case  and  the  water 
level  was  varied  in  the  other  case.  The  cycle  of  variation  was  20  time 
steps  for  both  of  the  hypothetical  simulations.  The  incident  wave  angle 
was  set  30  deg  in  the  most  seaward  calculation  cell  (end  of  the  grid)  and 
the  wave  period  to  8  sec.  In  the  first  case  the  wave  height  was  constant  1 
m  and  the  water  level  had  an  amplitude  of  0.5  m.  For  the  second  case  the 
water  level  was  fixed  but  the  wave  height  varied  between  0.5  and  1.5  m. 

108.  Figure  10a  displays  the  beach  profiles  at  the  last  time  step 
for  selected  profile  lines  along  the  grid  (profile  No.  1,  2,  3,  5,  and  10 
counting  in  the  direction  away  from  the  groin) .  Significant  accumulation 
occurred  along  the  profile  line  immediately  updrift  the  groin  whereas 
further  away  from  the  groin  the  accretion  was  less  marked.  The  steepening 
of  the  profile  in  the  seaward  part  closest  to  the  groin  is  in  qualitative 
agreement  with  what  is  observed  in  the  field.  Also,  note  that  the  bar 
feature  which  is  present  in  the  profiles  far  away  from  the  groin  is  almost 
completely  obliterated  in  the  profile  adjacent  to  the  groin. 

109.  In  Figure  10b  is  the  wave  height  and  longshore  current  dis¬ 
tribution  at  the  initial  and  final  time  step  for  the  profile  closest  to  the 


Profile  Depth  (m) 


0  50  100  150  200 

Cross-shore  distance  (m) 


-  Profile  1  -  Profile  2  Profile  3 

Profile  5  -  Profile  10 


(a)  beach  profiles  at  final  time  step  for  profile  lines  updrift  a  groin 
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(b)  beach  profile  ac  initial  and  final  time  step  and  corresponding 
wave  height  and  longshore  current  distribution  at 
profile  closest  to  groin 


Figure  10.  Beach  change  for  the  case  with  constant  wave  height  and  sinusoid 
ally  varying  water  level 
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groin  shown.  The  waves  shoal  up  and  break  close  Co  shore  at  the  initial 
profile  (the  break  point  located  at  about  60  m  from  an  arbitrary  base 
line),  whereas  the  break  point  moves  considerably  offshore  when  accumula¬ 
tion  takes  place  (break  point  at  about  125  m) .  The  longshore  current 
distribution  is  peaked  for  the  initial  profile  but  becomes  very  flat  as 
material  accretes  updrift  the  groin. 

110.  At  Figures  11a  and  lib  is  the  corresponding  beach  evolution 
displayed  for  the  case  of  a  constant  water  level  and  a  sinusoidally  varying 
wave  height.  The  beach  change  is  very  similar  as  compared  to  the  case  with 
varying  water  level  although  the  bar  feature  is  even  less  noticeable 
immediately  adjacent  to  the  groin.  The  accumulation  updrift  the  groin  is 
stronger  and  the  seaward  movement  of  the  break  point  is  larger.  Note  that 
wave  reformation  takes  place  at  the  last  time  step  and  two  breakpoints 
occur  in  Figure  lib. 
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(a)  beach  profiles  at  final  time  step  for  profile  lines  updrift  a  groin 
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(b)  beach  profile  at  initial  and  final  time  step  and  corresponding 
wave  height  and  longshore  current  distribution  at 
profile  closest  to  groin 


Figure  11.  Beach  change  for  the  case  with  constant  water  level  and  sinusoid 
ally  varying  wave  height 
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c*********************************************************************** 

C*  PROGRAM  3DBEACH 

C*********************************************************************** 
C*  IS  A  3D -MODEL  FOR  SIMULATING  SHORELINE  AND  BEACH  PROFILE  EVOLUTION 
C******************************* **************************************** 

c* 

C*  PROGRAM  EDITORS: 

C*  HANS  HANSON  MAGNUS  LARSON 

C*  DEPARTMENT  OF  WATER  RESOURCES  ENGINEERING 

C*  UNIVERSITY  OF  LUND,  SWEDEN 

C*  AND 

C*  NICHOLAS  C.  KRAUS 

C*  COASTAL  ENGINEERING  RESEARCH  CENTER 

C*  U.S.  ARMY  ENGINEER  WATERWAYS  EXPERIMENT  STATION 
C*  VICKSBURG,  MS,  USA 

C*  AUGUST  1988. 

C* 

C*********************************************************************** 

c* 

C***  THIS  PROGRAM  CALCULATES 

c***  -  THE  BREAKING  WAVE  HEIGHT  AND  ANGLE 
C*** 

c***  -  WAVE  HEIGHT  DECAY  INSIDE  BREAKER  ZONE 
C*** 

C***  -  CROSS -SHORE  TRANSPORT  RATES  QX(Y) 

C***  -  LONGSHORE  CURRENT  VELOCITY  DISTRIBUTION  V(X) 

C*** 

c***  -  LONGSHORE  TRANSPORT  RATE  DISTRIBUTION  QY(X) 

C*** 

C***  -  NEW  BOTTOM  BATHYMETRY  DEP(X.Y)  USING  QX(Y)  AND  QY(X) 

C*** 

C*** 

C***  INDEX:  Z-DEEP  WATER,  IN-INITIAL,  BR-BREAKING 

C***  SWL-SEAWALL,  GR-GROIN,  OB-OFFSHORE  BOUNDARY 

C*** 

c *********************************************************************** 

c*** 

C* . 

PARAMETER(MM-100 , NN-50 , NNGR-10 , NM-NN*MM , 

&M1-MM+1 , Nl-NN+1 , NM1-M1*N1 ) 

REAL  KX.KY.NY 
CHARACTER*70  TITLE 
C* 

REAL  D(MM,NN) , DSMO(MM.NN) , E(M1) ,QXT(M1 ,N1) 

REAL  H(M1) ,QX(M1,N1) ,QY(M1,N1) , V(M1) , DISS (MM) 

REAL  DIN(MM.NN) ,CEXO(5) ,ZC(MM) , QYX(M1 , Nl) , DCHG 
REAL  HBEG,ZBEG,T,DHTOT(Ml) ,DH(M1) ,DMEAN(MM) 

REAL  QSMOX(Ml) ,QSMOY(Ml) , VPRI (Ml) , HPRI (Ml ) 

INTEGER  ISHORE(NN) ,NBR(M1) , IYGR(NNGR) 

INTEGER  IGTIP(NNGR) , ISEAWL(NN) 


\ 
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c* 

COMMON  /BLOCKA/DX , DY , DTX , M , N , G , PI , GAMMA 
COMMON  /BL0CKB/D50 , SRATIO , NY , SLOPE , DFS 
COMMON  /BLOCKC/CQ , CEQ , CEXO , EPS 
COMMON  /BLOCKD/CRED , C  EXR , AEXP 
COMMON  /BLOCKE /HBEG , ZBEG.T, DCHG , CC , DOFF 
COMMON  /BLOCKF/CF,GAMMIX,KY,QCUT 
COMMON  /BLOCKG/IBR1GH, IBLEFT 
C* 

C*INITIALIZE 

c* 

DATA  QX/NM1*0 . / , QY/NM1*0 . / , V/M1*0 . / 

DATA  G/9. 806/, PI/3. 141593/, GAMMA/1.0/ 

DATA  NY/1 . 2E-6/ , SRATIO/2 . 65/, SLOPE/O . 02/ , EPS/2 . E- 3/ 

DATA  CRED/0 . 2/ , CEXR/0 . 1/ , AEXP/0 . 5/ , DFS/0 . 5/ 

DATA  QCUT/0 . 0E- 5/CC/1226 . / 

C* 

C*FILE  PREPARATION  6,  READING 
C* 

OPEN  (UNIT-10, NAME- 'START. DAT' , STATUS- ' UNKNOWN ' ) 

OPEN  (UNIT-12 .NAME- 'DEPTH. DAT' , STATUS- 'UNKNOWN ' ) 

OPEN  (UNIT-16, NAME-'WAVES. DAT’ , STATUS-' UNKNOWN • ) 

OPEN  (UNIT-1 8, NAME- 'TRANY.DAT' , STATUS- 'UNKNOWN ' ) 

OPEN  (UNIT-20, NAME- 'TRANX.DAT' , STATUS- 'UNKNOWN' ) 

OPEN  (UNIT-22, NAME-'DEPUT. DAT' , STATUS- 'UNKNOWN ' ) 

C* 

C*  READ  FILE  HEADING 
C* 

DO  8  1-1,4 
READ(10,*) 

8  CONTINUE 
C* 

C*  RUN  TITLE 
C* 

READ( 10 , *) 

READ( 10 , 503 )  TITLE 
C* 

C*  NUMBER  OF  CALCULATION  TIME  LOOPS  &  GRID  CELLS 
C* 

READ(10 ,*) 

READ( 10 ,*)  NCL.M.N 
C* 

C*  ERROR  CHECK  AND  MESSAGE 
C* 

IF(N.GT.NN)  THEN 

WRITE(*,*)  'ERROR.  TOO  MANY  CALCULATION  CELLS  IN  X-DIRECTION 
WRITE(* , *)  'PLEASE  CHANGE  VALUE  OF  "N"  IN  INPUT  FILE.’ 
NOCON-1 
GOTO  116 
END  IF 

IF(M.GT.MM)  THEN 

WRITE(* , *)  'ERROR.  TOO  MANY  CALCULATION  CELLS  IN  Y- DIRECTION 
WRITE(* , *)  ’PLEASE  CHANGE  VALUE  OF  "M"  IN  INPUT  FILE.’ 
NOCON-1 
GOTO  116 
END  IF 
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i 

\ 


) 

\ 
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c* 

C*  CROSS -SHORE  AND  LONGSHORE  TIME  INCREMENT  &  SPACE  INTERVAL 
C* 

READ ( 10 , *) 

READ ( 10 , *)  DTX , DX , DY 

DTXM-DTX 

DTX— DTX*60 . 

C* 

C*  LONGSHORE  TRANSPORT  PARAMETERS 
C* 

READ( 10 , *) 

READ (10 ,*)  CF.GAMMIX 
C* 

C*  MEDIAN  SAND  GRAIN  SIZE  (IN  MM) 

C* 

READ (10,*) 

READ (10,*)  D50 
D50-D50/1000. 

C* 

C*  OUTPUT  INTERVAL  EXPRESSED  AS  NO.  OF  CALCULATION  LOOPS 
C* 

READ ( 10 , *) 

READ ( 10 , *)  I OUT 
C* 

C*  DEPTH  CORRESPONDING  TO  OFFSHORE  WAVE  DATA 
C* 

READ (10 ,*) 

READ (10,*)  DOFF 
C* 

C*  SAND  TRANSPORT  CALIBRATION  COEFFICIENTS  (KY-L-S,  KX-C-S) 

C* 

READ(10,*) 

READ(10,*>  KY.KX 
CQ— KX* 1 . E  -  6 

c* 

C*  INPUT  BOUNDARY  CONDITIONS  AT  RIGHT  AND  LEFT  SIDE  OF  GRID 
C*  (O-OPEN  COAST,  1-IMPERMEABLE  GROIN) 

C* 

READ (10  ,*) 

READ ( 10 , *)  IBRIGH, IBLEFT 
C* 

C*  INITIAL  SCREEN  MESSAGES 
C* 

WRITE(* ,*)  CHAR(12) 

WRITE(*,*)  'COASTAL  ENGINEERING  RESE 
&H  CENTER’ 


WRITE(*,*) 
WRITE!*,*) 
WRITE(* , *) 

'  ***** 

****** 

****** 

******* 

***** 

** 

WRITE(*,*) 

** ' 

•  ******* 

******* 

******* 

******* 

******* 

&***  ** 
WRITE(* , *) 

** ' 

•  **  ** 

**  ** 

**  ** 

** 

**  ** 

&  **  ** 
WRITE!* , * ) 

**  ' 

'  **  ** 

**  ** 

**  ** 

** 

**  ** 

&  **  irk 

**  ' 

r 

i 


\ 


R  C 

*+* 

**** 

*> 
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WRITE(* . *) 

'  ** 

** 

** 

** 

** 

** 

** 

** 

** 

&  ** 

**' 

WRITE(* , *) 

•  ** 

** 

** 

****** 

***** 

******* 

** 

&  ******* ' 

WRITE(*,*) 

'  ** 

** 

** 

****** 

***** 

******* 

** 

&  ******* ’ 

WRITE(* ,*) 

'  ** 

** 

** 

** 

** 

** 

** 

** 

** 

&  ** 

**  * 

WRITE(* ,*) 

*  ** 

** 

** 

** 

** 

** 

** 

** 

** 

&  **  ** 

**' 

WRITE(*,*) 

'** 

** 

** 

** 

** 

** 

** 

** 

** 

&  **  ** 

**' 

WRITE(*,*) 

'  ******* 

******* 

******* 

******* 

** 

** 

**** 

&***  ** 

**' 

WRITE(*,*) 

•  ******* 

****** 

****** 

******* 

** 

** 

*** 

&**  ** 

**' 

WRITE(* ,*) 

C* 

C*  CONFIRMATION  ON  SCREEN 
C* 

WRITE(* , 612 )  'RUN:  ’.TITLE 
WRITE(*,*) 

VJRITE(* , *)  '  KX  ='  ,KX,  '  KY  — '  , KY , '  NCL  -' ,NCL 

WRITE(* , * )  '  DX-’.DX,'  DTX  - ' , DTXM , '  MIN' 

WRITE(* , *)  '  M  -' ,M, '  N  -'  ,N 
C* 

C*  READ  AND  SAVE  INITIAL  BOTTOM  CONFIGURATION 
C* 

CALL  DEPIN(DIN) 

DO  29  I-l.M 
DO  30  J-l.N 

D(I ,J)-DIN(I ,J) 

30  CONTINUE 

29  CONTINUE 
C* 

C*  CALCULATING  EQ.  PROFILE  SHAPE  FACTOR 
C* 

DD50-D50*1000. 

IF(DD50 . LT . 0 . A)  ADEAN-0 . 41*DD50** . 94 
I F (DD50 . GE . 0 . 4 . AND . DD50 . LT . 10 . )  ADEAN-0 . 23*DD50** . 3 2 
IF(DD50.GE. 10. . AND . DDSO . LT . 40 . )  ADEAN-0 . 23*DD50** . 28 
IF(DD50.GE.40. )  ADEAN-0 ,46*DD50** . 11 
C* 

C*  EQUILIBRIUM  ENERGY  DISSIPATION 
C* 

CEQ-0 . 75*ADEAN**1 . 5*3902 . 

C* 

C ***** ************************ 

C*  MAIN  CALCULATION  LOOP  * 

C***************************** 

C* 

DO  115  NS-l.NCL 
KKK-0 

IF(NS . EQ . 1 . OR . NS . EQ. NCL)  KKK-1 
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c* 

C*  OBTAIN  WAVE  DATA  ALONG  THE  COAST  AT  END  OF  GRID 
C* 

CALL  WAVIN(D) 

C* 

C*  ADD  WATER  LEVEL  CHANGE 
C* 

DO  49  J-l.N 
DO  48  I-l.M 

D(I,J)-D(I,J)+DCHC 

48  CONTINUE 

49  CONTINUE 
C* 

C*  SMOOTH  BOTTOM  PROFILE  FOR  WAVE  CALCULATION 
C* 

MAV-9 

DO  51  J-l.N 

DO  52  I-M ,1,-1 

IF(I . GT . M-MAV/2 )THEN 
DSMO(I,J)-D(I,J) 

ELSEIF ( I . LE .M-MAV/2 .AND . D( I -MAV/2 , J ) . GT . 0)THEN 
SUM=0 

DO  53  K-I -MAV/2 , I+MAV/2 
SUM-SUM+D ( K , J ) 

53  CONTINUE 

DSMO (I , J) -SUM/REAL(MAV) 

ELSE 

DO  54  K-1,I 

DSMO(K,J)-D(K,J) 

54  CONTINUE 
GOTO  51 

ENDIF 

52  CONTINUE 

51  CONTINUE 
C* 

C*  CALCULATE  FOR  EACH  PROFILE 

C*  _ 

C* 

DO  87  J-l.N+l 
C* 

C*  DETERMINE  MEAN  PROFILE 
C* 

DO  68  I-l.M 

IF(J.EQ.1)THEN 
DMEAN (I)-D(I,J) 

ELSEIF(J.EQ.N+I)THEN 
DMEAN (I)-D(I.J-l) 

ELSE 

DMEAN(I)— 0. 5*(D(I,J)+D(I,J-1)) 

ENDIF 

68  CONTINUE 

C* 

C*  DETERMINE  THE  WAVE  HEIGHT  DISTRIBUTION 
C* 

CALL  WAVEHD(KKK , DSMO , J , H , E, V , DISS , I  STOP , NBR , ZC , 
6.  ZCM) 


\ 


56 


I 


DO  69  I-l.M+l 

I F ( I . EQ . 1 ) THEN 
DH  ( I )  -DMEAN  ( I ) 

ELSEIF(I.EQ.M+1)THEN 
DH(I)— DMEAN(M) 

ELSE 

DH ( I ) — 0 . 5* ( DMEAN ( I ) +DMEAN ( I  - 1 ) ) 

END  IF 

DHTOT ( I ) — DH (I)+E(I) 

69  CONTINUE 

C* 

C*  DETERMINE  CROSS -SHORE  TRANSPORT 
C* 

CALL  TRANX ( J , H , DISS , NBR , DMEAN , DH , DHTOT , NRU , 

&  NFS ,QXT) 

C* 

C*  DETERMINE  LONGSHORE  TRANSPORT 
C* 

CALL  TRANY ( J , I STOP , NRU , NFS , NBR , V , H , ZC , ZCM , QY , QYX ) 

C* 

C*  WRITE  WAVE  HEIGHT  AND  LONGSHORE  CURRENT  AT  FIRST  AND  LAST  TIME 
C*  STEP  FOR  PLOTTING 
C* 

IF(KKK. EQ. 1 . AND. J . EQ. 1)THEN 
NPRI=ISTOP 
DO  383  I-NPRI.M 
VPRI(I)=-V(I) 

HPRI ( I)=H( I ) 

383  CONTINUE 

HPRI (M+1)-H(M+1) 

END  IF 

87  CONTINUE 

C* 

C*  SPECIAL  CASE  OF  OPEN  COAST  AT  RIGHT  END  OF  GRID 
C* 

IF(IBRIGH . EQ . 0)THEN 
DO  38  1=1, M 

QY ( I , 1 ) =QY (1,2) 

38  CONTINUE 

END  IF 
C* 

C*  INTERPOLATE  CROSS -SHORE  TRANSPORT  RATES  ALONG  DEPTH  GRID 
C* 

DO  33  J-l.N 

QX(1 , J )=0 . 0 
DO  43  1=2. M 

QAA-0 . 2  3* (QYX ( I , J ) +QYX ( I  - 1 , J ) +QYX ( I , J+l ) 

&  +QYX(I-1,J+1)) 

QX(I ,J)=0.5*(QXT(I,J)+QXT(I ,J+1))+QAA 
43  CONTINUE 

QX(M+1,J)=0.0 
33  CONTINUE 
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c* 

c*  change  back  to  original  water  level  before  applying  the 

C*  MASS  CONSERVATION  EQUATION  . . 

C* 

DO  57  J-1,N 
DO  58  I-l.M 

D(I,J)-D(I,J)-DCHG 

C* 

C*  SMOOTH  TRANSPORT  RATES  ALONGSHORE  AND  CROSS -SHORE 
C* 

IF( I . GT. 1 .AND , I . LT .M)THEN 

QSMOX(I )-0 . 25*QX(I - 1 , J )+0 . 50*QX(I , J )+0 .  25*QX( 1+1  J) 
QSMOY(I )-0 . 25*QY ( I  - 1 , J )+0 . 50*QY(I , J)+0 . 25*QY( I+I  J) 
ENDIF 

58  CONTINUE 

DO  100  I-l.M 

QX ( I , J ) -QSMOX ( I ) 

QY ( I , J ) -QSMOY ( I ) 

100  CONTINUE 

57  CONTINUE 

C* 

C*  DETERMINE  THE  CHANGES  IN  BATHYMETRI 
C* 

CALL  BATHYM(QX,QY,D) 

IF(KKK. EQ. 1)THEN 
WRITE( 50 , *)  M 
DO  333  I-l.M 

WRITE( 50 , *)  REAL( I - 1)*DX, -D( I  1) 

CONTINUE 

WRITE(50.*)  M-NPRI+1 
DO  334  I-NPRI.M 

WRITE(50 ,*)  REAL(I-1)*DX  VPRI(I) 

CONTINUE 

WRITE(50,*)  M-NPRI+1 
DO  335  I-NPRI.M 

WRITE(50,*)  REAL(I-1)*DX.0.5*(HPRI(I)+HPRI(I  +  1)  ) 
CONTINUE 
DO  336  1=1, M 

WRITE(51 , 110)  - D ( I , 1) , VPRI (I) ,0 . 5*(HPRI (I)+HPRI (1+1)1 
FORMAT ( '  ' , 3F8 . 3) 

CONTINUE 
ENDIF 
C* 

C*  INTERMEDIATE  OUTPUT 
C* 

IF(MOD(NS , IOUT) , EQ . 0 . OR, NS . EQ.NCL)  THEN 

c* 

C*  WRITE  L~S  TRANSPORT  QY,  C-S  TRANSPORT  QX,  AND  DEPTH 

c* 

IMAT-N/10 

IF(IMAT.NE.O)  THEN 
IMOD-MOD(N.IO) 

II-O 

IIM--9 

DO  50  IMOUT-1 , IMAT 
II-II+10 


333 


334 


335 


110 

336 
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IIM-IIM+10 

WRITE(18 ,*)  'QY*10**4  FOR  J-  TO  ',11, 

&  '  FOR  TIMESTEP  ' ,NS 

WRITE(18, 600)  ((QY(I,J)*10. **4 , J-IIM , I I ) , 1-1 ,M) 
WRITE(20 , *)  ’QX*10**4  FOR  J-  TO  ’.II, 

&  '  FOR  TIMESTEP  ' ,NS 

WRIT£(20, 600)  ( (QX(I ,J)*10.**4, J-IIM, II) , 1-1 ,M) 

WRITE(22,*)  'DEPTH  FOR  J-  ',IIM,'  TO  ',11, 

&  '  FOR  TIMESTEP  ' ,NS 

WRITE(22,600)  ( (-D(I.J), J-IIM, II), I-l.M) 

50  CONTINUE 

ENDIF 

IF(IMOD.NE.O)  THEN 
II=II+IMOD 
IIM-IIM+10 
WRITE(18,*) 

WRITE(18 ,*)  'QY*10**4  FOR  J-  ',IIM,'  TO  '.II 
WRITE( 18 , 600)  ( (QY(I,J)*10.**4, J-IIM, II) ,1-1, M) 
WRITE(20 , *)  ’QX*10**4  FOR  J-  ',IIM,'  TO  ',11 
VJRITE( 20 , 600)  ( ( QX ( I ,  J )*10 .**4 ,  J-IIM,  II)  ,  1-1  ,M) 
WRITE(22,*)  'DEPTH  FOR  J-  ',IIM,'  TO  ',11 
WRITE(22 , 600)  ( (-D(I,J), J-IIM, II), I-l.M) 

ENDIF 

ENDIF 

115  CONTINUE 

116  IF(NOCON.EQ.l)  THEN 

WRITE(* , *)  'EXECUTION  STOPPED' 

ENDIF 

C* 

C*  END -OF -COMPUTATION  SIGNAL  (BELL) 

C* 

WRITE(* , *)  CHAR ( 7 ) , CHAR ( 7 ) 

C* 

C*  FORMAT  STATEMENTS  FOR  OUTPUT 
C* 

503  FORMAT (3X.A70) 

600  FORMAT ( 10F7 . 3 ) 

601  F0RMAT(/1X, 'OUTPUT  LAST  TIME  STEP  NO.', 17/) 

602  FORMAT (IX, 1017) 

603  FORMATC/1X, 'BREAKING  WAVE  HEIGHT'/(1X, 10F7.2)) 

604  F0RMAT(/1X, 'BREAKING  WAVE  ANGLE '/( IX , 10F7 . 1 ) ) 

605  F0RMAT(/1X, 'BREAKING  WAVE  LOCATION '/( IX , 1017 ) ) 

610  F0RMAT(10I7) 

612  FORMAT ( IX, A4, IX, A70) 

614  FORMAT (1X,A,I3,A,I3) 

616  FORMAT (10F7. 2) 

END 
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c* 

c********************************************************************** 

SUBROUTINE  DEPIN(DEP) 

C  *★************★*****★*********★★*★**★**********■***■*'*  ■*■*****■***•*  •*•*■***«** 

c* 

C*  READS  DEPTHS  FROM  A  FILE 
C* 

PARAMETER(MM-100 , NN-50 , NNGR-10 , NM-NN*MM, 

&M1-MM+1 , Nl-NN+1 , NM1-M1*N1 ) 

C* 

REAL  DEP(MM.NN) 

C* 

COMMON  /BLOCKA/DX , DY , DTX , M , N , G , PI , GAMMA 
C* 

C*  READ  DEPTHS  ALONG  THE  COAST 
C* 

IMAT-N/10 

IF(IMAT.NE.O)  THEN 
IMOD— MOD(N ,  10) 

II-O 

IIM--9 

DO  5  IDEP-l.IMAT 
II-II+10 
IIM-IIM+10 
READ(12 ,*) 

READ(12, 600, END-99)  ((DEP(I.J) ,J-IIM,II) ,1-1, M) 

5  CONTINUE 

END  IF 

IF(IMOD.NE.O)  THEN 
II-II+IMOD 
IIM-IIM+10 
READ ( 12 , *) 

READC12, 600, END-99)  ( <DEP(I , J) , J-IIM, II)  ,  1-1  ,M) 

ENDIF 
GOTO  60 

99  CONTINUE 

WRITE(*,*)  'ERROR  FOUND  IN  "DEPIN" . ' 

WRITE (*,*)  'FILES  "DEPTH"  (AND  "WAVES")' 

WRITE(*,*)  'CONTAIN  TOO  FEW  VALUES.' 

WRITE(* , *)  'PLEASE  CHANGE  IN  DATA  FILE(S) . ' 

NOCON-1 
GOTO  100 
60  CONTINUE 

C* 

600  FORMAT (10F7. 2) 

100  CONTINUE 
RETURN 
END 
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c* 

c**************************************************************  •*****■**■*■*■ 
SUBROUTINE  WAVIN(D) 

C *********************************************************************** 

c* 

C*  READS  OFF-SHORE  WAVE  DATA  FROM  A  WAVE  FILE 
C* 

PARAMETER (MM-lOO , NN-50 , NNGR-10 , NM-NN*MM , 

&M1-MM+1 , Nl-NN+1 , NM1-M1*N1) 

REAL  D(MM,NN) ,HBEG, ZBEG, DCHG 
COMMON/BLOCKA/DX , DY , DTX , M , N , G , PI , GAMMA 
COMMON/BLOCKE/HBEG , ZBEG , T , DCHG , CC , DOFF 
C* 

C*  READ  DEEP  WATER  WAVE  CONDITIONS 
C* 

66  READ(16,*, END-99, ERR-77)  HOFF, ZOFF.T, DCHG 

HBEG-HOFF 
ZBEG— ZOFF 
C* 

C*  REACHED  END  OF  WAVE  DATA  FILE.  REWIND,  READ  HEADING  AND  CONTINUE. 

C* 

RETURN 

99  REWIND  16 

GOTO  66 

77  WRITE(*,*)  'ERROR  IN  WAVIN’ 

STOP 

END 


I 


c* 

C**********************^************************************'*’*“Ar****’*,**'A-'5!f 

SUBROUTINE  TRANX ( IP, H , DISS , NBR , DMEAN , DH , DHTOT , NRU , NFS , QX) 
C*********************************************************************** 

c* 

C*  THIS  ROUTINE  CALLS  THE  SUBROUTINES  NECESSARY  FOR  GENERATING 
C*  THE  CROSS-SHORE  TRANSPORT  RATES  FOR  A  SPECIFIC  PROFILE  J 
C* 

PARAMETER (MM-100 , NN-50 , NNGR-10 , NM-NN*MM, 

&M1-MM+1 , Nl-NN+1 , NM1-M1*N1) 

REAL  QX(M1,N1) , DX, DY, DTX, Q(M1) , DMEAN (MM) , 

&H(M1) ,DISS(MM) ,DH(M1) .DHTOT(Ml) 

INTEGER  NBR (Ml) ,NPE(5) , NPP(5) , NPB(5) , SIGN ,NRU , NFS 
C* 

COMMON  /BLOCKA/DX , DY , DTX , M , N , G , PI , GAMMA 
C* 

C*  DETERMINE  LOCATIONS  OF  TRANSPORT  REGIONS 
C* 

CALL  DOMAIN ( DH , DHTOT , NBR , H , NRU , NFS , NPE , NPP , NPB , NRBPE) 

C* 

C*  DETERMINE  DIRECTION  OF  TRANSPORT 
C* 

CALL  TRSIGN(SIGN) 

C* 

C*  DETERMINE  CROSS -SHORE  TRANSPORT  RATES 
C* 

CALL  XTRANS ( SIGN , DMEAN , NRU , NFS , NPE , NPP , NPB , NRBPE , DISS , Q) 

DO  20  I-l.M+l 
QX(I , IP)-Q(I) 

20  CONTINUE 
RETURN 
END 
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c* 

c*********************************************************************** 
SUBROUTINE  DOMAIN (DH , DHTOT , NBR , H , NRU , NFS , NPE , NPP , NPB , NRBPE) 

C****************************************************^**********'>'*-**^*** 

C* 

C*  THIS  ROUTINE  DETERMINES  THE  LOCATION  OF  THE  BOUNDARIES  OF 
C*  THE  DIFFERENT  CALCULATION  DOMAINS. 

C* 

INTEGER  MM , NN , M , N , I , NFS , NPP (5) ,NRU,NPE(5) ,NPB(5) 

PARAMETER (MM- 100 , NN-50 , NNGR-10 , NM-NN*MM , 

&M1-MM+1 ,  Nl-NN+1 ,  NM1-M1*N1) 

INTEGER  NBR(Ml) ,J,L, NRBPE 

REAL  DHTOT (MI) , DFS , LO , T , HO , DX, DY, DTX , EPS , GAMMA, DH (Ml) , 

&H (Ml ) , CQ , CEQ , DRU , SRATIO , D50 , CEXO ( 5 ) , NY , DCHG 
LOGICAL  0K1.0K2 
C* 

COMMON/BLOCKA/DX , DY , DTX , M , N , G , PI , GAMMA 
COMMON/BLOCKB/D50 , SRATIO , NY , SLOPE , DFS 
COMMON/BLOCKC/CQ , CEQ , CEXO , EPS 
COMMON /BLOCKE/HBEG , ZBEG , T , DCHG , CC , DOFF 
C* 

HO-HBEG 
DO  3  1-1,5 
NPP(I)-0 
NPE(I)-0 
NPB(I)-0 
3  CONTINUE 

LO-1 . 5613*T**2 

J-0 

L-0 

DO  5  I-l.K-l 
C* 

C*  DETERMINE  LOCATION  OF  BREAK  POINTS,  PLUNGE  POINTS  AND 
C*  CORRESPONDING  SPATIAL  DECAY  COEFFICIENTS 
C* 

IF(NBR( I )  EQ. 1 . AND.NBR(I+1) . EQ . 0)THEN 
J-J+l 
NPB(J)—I 

CEXO( J)=0 . 40*(D50*1E3/H( I ) )**0 .47 
NPP( J)=NPB(J ) -NINT( 3*H ( I )/DX) 

ENDIF 

C* 

C*  DETERMINE  LOCATION  OF  REFORMATION  POINTS 
C* 

IF(NBR( I) . EQ . 0 . AND ,NBR(I+1) . EQ . 1)THEN 
L— L+l 

NPE(L)=I+1 
ENDIF 
5  CONTINUE 
NRBPE-J 
C* 

C*  DETERMINE  THE  RUNUP  LIMIT  ACCORDING  TO  RELATIONSHIP  FROM  CE 
C*  AND  CRIEPI  DATA 
C* 

DRU--1 ,47*HO*(TAN(SLOPE)/SQRT(HO/LO))**0. 79 
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c* 

C*  DETERMINE  THE  LOCATION  OF  THE  SURF  ZONE  END  AND  RUNUP  LIMIT. 
C*  WHEN  A  BOUNDARY  DEPTH  OCCUR  IN  THE  MIDDLE  OF  A  GRID  CELL  THE 
C*  THE  SEAWARD  BOUNDARY  NUMBER  OF  THE  CELL  IS  ASSIGNED  TO  THE 
C*  BOUNDARY  INTEGER  NUMBER 
C* 

OKI-. TRUE. 

OK2- . TRUE . 

DO  9  I— M ,1,-1 

IF((DH(I) .LT. DRU. AND . DH(I+1) .GE. DRU) . AND . OKI ) THEN 
NRU-I+1 
OKI-. FALSE. 

END  IF 

IF((DHTOT(I) . LT . DFS . AND. DHTOT(I+l) . GE . DFS) . AND. OK2) THEN 
NFS-I+1 
OK2-. FALSE. 

ENDIF 

IF( ( . NOT . OKI) .AND . ( . NOT . OK2) )  GOTO  88 
9  CONTINUE 

IF(OKl)  NRU=1 
88  DO  20  1-1 , NRBPE 

IF(NPEd) . GT . NPP( I ) )THEN 
NPP(I)-NPE(I) 

ENDIF 
20  CONTINUE 

RETURN 
END 
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c* 

c*********************************************************************** 

SUBROUTINE  TRSIGN(SIGN) 

C*********************************************************************** 

c* 

C*  THIS  ROUTINE  DETERMINES  THE  DIRECTION  OF  THE  SAND 
C*  TRANSPORT  ALONG  THE  PROFILE.  THE  TRANSPORT  IS  DIRECTED 
C*  ONSHORE  OR  OFFSHORE  DEPENDING  ON  THE  VARIABLE  SIGN 
C*  (SIGN  -  X,  OFFSHORE,  -  -1,  ONSHORE)  ACCORDING  TO  AN 
C*  EMPIRICAL  EXPRESSION. 

C* 

INTEGER  SIGN 

REAL  T , D50 , LO , HO , VF , WSTP , DIMVF , SRATIO , B , NY , DCHG 
COMMON/BLOCKB/D50 , SRATIO , NY , SLOPE , DFS 
COMMON/BLOCKE/HBEG , ZBEG , T , DCHG , CC , DOFF 
HO-HBEG 
LO-1 . 56I3*T**2 
WSTP-HO/LO/1 . 6 
C* 

C*  DETERMINE  THE  FALL  VELOCITY  ACCORDING  TO  HALLERMEIER 
C* 

B=(SRATIO- 1 )*9 . 81*D50**3/NY**2 
IF(B . LE . 39 ) THEN 

VF-(SRATIO-l)*9 . 81*D50**2/18/NY 
ELSEIF(B . GT . 39 . AND . B . LE . 1E4)THEN 

VF-( ( SRATIO- 1 )*9 . 81)**0 . 7*D50**1 . l/6/NY**0 . 4 
ELSE 

VF-( (SRATIO- 1)*9 .81*D50/0 . 91)**0 . 5 
END  IF 

DIMVF-HO/T/VF/I . 6 
C* 

C*  USE  EMPIRICALLY  DERIVED  CRITERIA  TO  PREDICT  TRANSPORT  DIRECTION 
C* 

IF(WSTP.LT.7E-4*DIMVF**3)THEN 

SIGN-1 

ELSE 

SIGN—- 1 
ENDIF 
RETURN 
END 


65 


c* 

c**********************************************************************-* 

SUBROUTINE  XTRANS (SIGN, DMEAN , NRU , NFS , NPE , NPP , NPB , NRBPE , DISS , Q ) 

C  ★*★*★*****★**★★**★**★★******★******★*★***★*★*★**★**★★■***★'********■*■*'*■*•** 

c* 

C*  THIS  ROUTINE  GENERATES  THE  CROSS -SHORE  TRANSPORT  RATES. 

C* 

INTEGER  I ,NFS ,NPP(5) ,NPE(5) ,  NZ  ,MM  ,NN  ,NRU  ,M,  SIGN ,  NPB(5 ) , 

&NRBPE.N 

PARAMETER (MM-100 , NN-50 , NNGR-10 , NM-NN*MM , 

&M1-MM+1 , Nl-NN+1 , NM1-M1*N1 ) 

REAL  DX.DTX, DMEAN (MM) ,Q(M1) ,QSMO(Ml) ,DISS(MM) ,CEXO(5) , DY , 

&CQ , CEQ , EPS , CRED , CEXR , AEXP , CEXP 
LOGICAL  END, REDUCE 
C* 

COMMON /BLOCKA/DX , DY, DTX,M,N,G , PI , GAMMA 
COMMON/BLOCKC/CQ , CEQ , CEXO , EPS 
COMMON/BLOCKD/CRED , CEXR , AEXP 
C* 

C*  DETERMINE  THE  Q- VALUES  ON  THE  NEXT  TIME  LEVEL 
C* 

END-. FALSE. 

IF(NPE(1) . LT.NFS)  NPE(1)=NFS 
IF(NPP(1) .LT.NFS)  NPP (1)— NFS 
C* 

C*  DETERMINE  TRANSPORT  RATES  IN  AREAS  OF  BROKEN  WAVES 
C* 

DO  5  I -1, NRBPE 
REDUCE-. FALSE. 

C* 

C*  IF  THE  BROKEN  WAVE  ZONE  IS  ONLY  ONE  CELL  WIDE  REDUCE  THE 
C*  TRANSPORT 
C* 

IF( I . EQ . 1 .AND . NPP ( I ) -NPE( I ) . LE . 1)  REDUCE- . TRUE . 

77  DO  10  J-NPE(I) ,NPP(I) 

IF( J . NE . NPE( I ) )THEN 

Q(J)»REAL(SIGN)*CQ*(0.5*(DISS(J)+DISS(J-1))-CEQ)+ 

&  EPS*(DMEAN( J ) - DMEAN (J -1 ) ) /DX 

IF(REAL(SIGN)*Q(J) . LT.0)THEN 
Q(J)-0 
END  IF 
ELSE 

IF(I . EQ . 1)THEN 

Q(J)-REAL(SICN)*CQ*(DISS(J)-CEQ)+ 

&  EPS* (DMEAN ( J+l) - DMEAN ( J ) ) /DX 

IF(REAL(SIGN)*Q(J) . LT.0)THEN 
Q(J)-0 
END  IF 
C* 

C*  REDUCTION  OF  TRANSPORT  RATE  FOR  NARROW  SURF  ZONE 
C* 

IF(REDUCE)  Q(J)-Q(J)*CRED 
ELSE 

Q(J)-REAL(SIGN)*CQ*(DISS(J) -CEQ)+ 

6.  EPS* (DMEAN  (  J  )  -  DMEAN (  J  - 1 ) ) /DX 

IF(REAL(SIGN)*Q(J) ,LT.0)THEN 
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Q(J)-0 
ENDIF 
END  IF 
ENDIF 

10  CONTINUE 
5  CONTINUE 

DO  25  I-l.NRBPE 

IF(I.EQ.NRBPE)  END-. TRUE. 

CEXP-CEXO(I)/5 . 

C* 

C*  DETERMINE  TRANSPORT  RATES  SEAWARD  OF  FIRST  BREAK  POINT 
C* 

IF(END)THEN 

C* 

C*  . . .  FROM  PLUNGE  POINT  TO  BREAK  POINT 
C* 

DO  11  J-NPP(I )+l , NPB( I ) 

Q(J)-Q(NPP(I))*EXP(-CEXP*DX*REAL(J-NPP(I) ) ) 

11  CONTINUE 
C* 

C*  .  .  .  FROM  BREAK  POINT  TO  END  OF  GRID 
C* 

DO  12  J-NPB(I)+1,M 

Q(J)=Q(NPB(I))*EXP(-CEXO(I)*DX*REAL(J-NPB(I)) ) 
IF(ABS(Q(J) ) . LT. ABS(Q(NPB(I) ) )/100 . )THEN 
NZ-J 
GOTO  66 
ENDIF 

12  CONTINUE 
ELSE 

C* 

C*  DETERMINE  TRANSPORT  RATES  BETWEEN  ZONES  OF  BROKEN  WAVES  SHOREWARD 
C*  OF  THE  EIRST  BREAK  POINT 
C*  ...FROM  PLUNGE  POINT  TO  BREAK  POINT 
C* 

DO  15  J=NPP( I )+l , NPB ( I ) 

Q(J)-Q(NPP(I) )*EXP( - CEXP*DX*REAL( J -NPP ( I ) ) ) 

15  CONTINUE 

C* 

C*  . . . FROM  BREAK  POINT  TO  REFORMATION  POINT 
C* 

NMP-NPB(I)+(NPE(I+1) -NPB(I) )/3 
C  WRITE(*,*)  NPP(I) ,NPB(I) ,NPE(I+1) ,NMP 

C  WRITE(*,*)  Q(NPP( I )) ,Q(NPB(I)) , Q(NPE( 1+1 ) ) 

DO  20  J=NPE( 1+1 ) - 1 , NMP , - 1 

Q( J )-Q(NPE( 1+1) )*EXP( -CEXR*DX*REAL(NPE( 1+1 ) - J ) ) 

20  CONTINUE 

DO  30  J-NPB ( I )  + 1 , NMP - 1 

Q(J)-(Q(NPB(I)) -Q(NMP) )*( 1- (REAL( J -NPB(I) )/ 

&  REAL (NMP- NPB (I ) ) )**AEXP)+Q(NMP) 

30  CONTINUE 

ENDIF 

66  END-. FALSE. 

25  CONTINUE 

IF(NPE(1) . NE . NFS)THEN 

DO  40  J-NPE( 1 ) - 1 , NFS , - 1 
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Q ( J ) — Q ( NPE  < 1 ) )*EXP( -CEXR*DX*REAL(NPE(1) -J) ) 
40  CONTINUE 

ENDIF 
C* 

C*  DETERMINE  TRANSPORT  RATES  IN  THE  SWASH  ZONE 
C* 

DO  50  J-NRU+1 , NFS  - 1 

Q(J)-Q(NFS)*(REAL( J ) -REAL(NRU) ) /REAL(NFS -NRU) 
50  CONTINUE 
C* 

C*  DETERMINE  TRANSPORT  RATES  AT  THE  BOUNDARIES 
C* 

Q(NRU)-0 

Q(NZ+l)-0 

C* 

C*  ZERO  Q- ARRAY  OUTSIDE  TRANSPORT  REGIONS 
C* 

DO  80  I-l.NRU-l 
Q(I)-0 

80  CONTINUE 

DO  90  I-NZ+2 , M 
Q(I)-0 

90  CONTINUE 

RETURN 
END 
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C* 

c**********************************************************-^********** 
SUBROUTINE  BATHYM(QX , QY , D) 

£***********************************************•******************•***** 

c* 

C*  THIS  ROUTINE  CALCULATES  THE  CHANGES  IN  BATHYMETRY  FROM 
C*  THE  MASS  CONSERVATION  EQUATION.  THE  TRANSPORT  RATE  MUST 
C*  BE  GENERATED  WHEN  THIS  SUBROUTINE  IS  CALLED 
C* 

INTEGER  MM,NN,I,J,M,N 

PARAMETER (MM=100 ,NN-50 .NNGR-XO ,NM-NN*MM, 

&M1-MM+1 .Nl-NN+l ,NM1-M1*N1) 

REAL  QX(M1,N1) ,D(MM,NN) , DX , DY , DTX ,QY(M1 ,N1) 

INTEGER  ISHORE(NN) 

C* 

COMMON /BLOCKA/DX , DY , DTX , M , N , G , PI , GAMMA 
C* 

DO  10  J-l.N 
DO  20  1=1, M 
C* 

C*  CALCULATE  THE  NEW  DEPTHS  ALONG  THE  GRID 
C* 

D(I,J)=D(I , J)+DTX/DX*(QX(I+1 , J) -QX(I,J))+ 

&  DTX/DY*(QY ( I , J+l) -QY( I , J) ) 

20  CONTINUE 

CALL  MAXANG ( J , D , 2 , M , I  SHORE ( J ) ) 

10  CONTINUE 
RETURN 
END 
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c* 

c*********-*-**************************-*********************************** 

SUBROUTINE  MAXANG( IP , D , NBEG.NEND , ISHOR) 
C******************************************************* **************** 

c* 

C*  THIS  ROUTINE  LIMITS  THE  STEEPNESS  OF  THE  PROFILE  BY  INITIATING 
C*  AVALANCHING  WHEN  THE  SLOPE  LOCALLY  EXCEEDS  A  PREDEFINED 
C*  MAXIMUM  VALUE.  THE  AVALANCHING  IS  CONSIDERED  TO  OCCUR  INFINITELY 
C*  FAST  IN  COMPARISON  WITH  THE  TIME  STEP  OF  THE  MODEL.  SAND  IS 
C*  REDISTRIBUTED  IN  THE  NEIGHBORING  CELLS  AT  A  SLOPE  CORRESPONDING 
C*  TO  THE  STABLE  SLOPE  AFTER  AVALANCHING. 

C* 

INTEGER  NBEG, NEND.MM.NN, I , NCORR.NDIR.NVAL, J , L, IMAX, IP, ISHOR , 

AM.N 

PARAMETER (MM=100 .NN-50 ,NNGR-10,NM-NN*MM, 

&M1-MM+1 ,  Nl-NN+1 ,  NM1— M1*N1 ) 

REAL  D(MM,NN) ,DD(MM) , BMAX , BAV , DX , DTX , BPMAX , DCORR, DCR , DSUM , 

&DDCRIT, DY 
LOGICAL  CHECK 

c* 

COMMON/BLOCKA/DX ,  DY ,  DTX , M , N , G ,  PI ,  GAMMA 
C* 

BAV-0 . 3 
BMAX-0 . 5 
L-0 
C* 

C*  DETERMINE  THE  MAXIMUM  SLOPE  ALONG  THE  PROFILE 
C* 

88  CHECK=. FALSE. 

L-L+l 

BPMAX-0 

DO  10  I=NBEG,NEND-1 

IF(ABS ( (D(X .IP) -D(I-l.IP) )/DX) . GT ,ABS( BPMAX) )THEN 
BPMAX-(D(I,IP)-D(I-1,IP))/DX 
IMAX-I 
ENDIF 
C* 

C*  DETERMINE  THE  LOCATION  OF  THE  SHORELINE  FOR  PROFILE  IP 
C* 

IF(D(I-1,IP) . LE . 0 . AND . D ( I , I P ) . GT . 0 )  ISHOR-I-1 
10  CONTINUE 
C* 

C*  CHECK  IF  THE  MAXIMUM  SLOPE  EXCEEDS  THE  ANGLE  OF  REPOSE 
C* 

IF (ATAN (ABS (BPMAX) ) . GT . BMAX)  CHECK-. TRUE. 

IF(CHECK)THEN 

WRITE(*  , *)  '/  MLANCHIIIIIII . ' 

C* 

C*  REDISTRIBUTE  THE  SAND  IN  THE  NEIGHBORING  CELLS  ACCORDING  TO  THE 
C*  SLOPE  AFTER  AVALANCHING  HAS  OCCURRED  (BAV) 

C* 

DCR*  DX*TAN ( BAV ) 

I  F(  BPMAX .  GT .  0 )  THEN 
NDIR--1 
ELSE 

NDIR-1 
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IMAX— IMAX-1 
END  IF 

NCORR-IMAX 

77  NCORR-NCORR+NDIR 

C* 

C*  CONTINUE  TO  REDISTRIBUTE  SAND  UNTIL  THE  FIRST  CELL  OUTSIDE  THE  AREA 
C*  WHERE  SAND  IS  MOVED  HAS  A  SLOPE  LOWER  THAN  BAV 
C* 

NVAL-1 

DSUM-0 

DO  20  I— IMAX+NDIR , NCORR , NDIR 

IF(I . GT . NEND . OR. I . LT . NBEG- 1)THEN 

WRITE(* ,*)  'AV  STOP I , IMAX+NDIR , NCORR , L 
STOP 
END  IF 

DSUM-DSUM+D( I , IP) 

NVAL=NVAL+1 
20  CONTINUE 

DCORR- - D ( IMAX ,  I P) * (NVAL- 1 ) /NVAL+DSUM/NVAL+DCR/2* ( NVAL- 1 ) 

DD ( IMAX) -D ( IMAX , I P) +DCORR 
J=1 
C* 

C*  CALCULATE  THE  NEW  DEPTHS 
C* 

DO  30  I-IMAX+NDIR, NCORR, NDIR 
J-J+l 

DCORR- DD( I MAX) -D(I , IP) -REAL(J- 1)*DCR 
DD(I)— D(I,IP) +DCORR 
30  CONTINUE 

I F (NCORR+NDIR . GE . NBEG - 1 ) THEN 
IF(NDIR.EQ. -1)THEN 

DDCRIT-(DD (NCORR) -D (NCORR+NDIR, IP) )/DX 
IF(DDCRIT.GT. TAN(BAV) )  GOTO  77 
ELSE 

DDCR I T- ( D ( NCORR+NDI R , I P ) • DD ( NCORR ) ) /DX 
IF(DDCRIT . LT, -TAN (BAV) )  GOTO  77 
ENDIF 
END  IF 

DO  40  I-IMAX, NCORR, NDIR 
D(I , IP)-DD(I) 

40  CONTINUE 

ENDIF 
C* 

C*  GO  BACK  AND  CHECK  IF  THE  ANGLE  OF  REPOSE  IS  EXCEEDED  IN  OTHER 
C*  PARTS  OF  THE  PROFILE 
C* 

IF(CHECK , AND , L. LT . 20)  GOTO  88 

IF(L.GE. 20)  WRITE(*,*)  'AVALANCHING  ROUTINE  NOT  CONVERGING’ 

IF(L. GE . 20)  STOP 

RETURN 

END 
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c************************************************************ 
SUBROUTINE  CONPLA ( D , J P , DMEAN , ZC , ZCM ) 
C************************************************************ 

c* 

C*  THIS  ROUTINE  CALCULATES  THE  AVERAGE  PROFILE  DEPTH  BETWEEN 
C*  TWO  NEIGHBORING  LINES.  ALSO,  THE  CONTOUR  ORIENTATION  IS 
C*  DETERMINED  USING  AN  AVERAGE  OF  TWO  PLANES  FITTED  THROUGH 
C*  THREE  POINTS  AROUND  THE  CALCULATION  POINT. 

C* 

PARAMETER  (MM-100 , NN-50 , NNGR-10 , NM-NN*MM , 

&M1-MM+1 , Nl-NN+1 , NM1-M1*N1 ) 

INTEGER  JL, JR, JP.KK 
REAL  D(MM ,NN) , DMEAN (Ml) , ZC(MM) , ZCM 
COMMON/BLOCKA/DX ,  DY ,  DTX ,  M ,  N ,  G ,  PI ,  GAMMA 
C* 

C*  CALCULATE  MEAN  WATER  DEPTH 
C* 

IF(JP . EQ . 1 ) THEN 
DO  10  I-l.M 

DMEAN (I)-D(I,1) 

10  CONTINUE 

JL=1 
JR-2 

ELSEIF(JP. EQ.N+1)THEN 
DO  20  I— 1 , M 

DMEAN ( I ) =D ( I , N ) 

20  CONTINUE 

JL-N-1 
JR-N 
ELSE 

DO  30  I-l.M 

DMEAN(I)-0.5*(D(I,JP)+D(I,JP-1)) 

30  CONTINUE 

JL-JP-1 
JR-JP 
END  IF 
C* 

C*  DETERMINE  CONTOUR  ORIENTATION 
C* 

DO  40  1-2, M-l 

IF(D(I+1,JL)-D(I-1,JL) .GT. 0.001) THEN 

BSC1-DX*(2*D(I,JR)-D(I-I,JL)-D(I+1,JL))/DY/(D(I+1,JL) 
&  -D(I-l.JL)) 

ELSE 

BSC1-1000. 

ENDIF 

Zl-ATAN(-BSCl) 

IF(D(I+1,JR)-D(I-1,JR) .GT. 0.001) THEN 

BSC2-DX*(D( 1  - 1 , JR)+D( 1+1 , JR) - 2*D( I , JL) )/DY/ (D(I  +  1,JR) 
&  -D( I  - 1 , JR) ) 

ELSE 

BSC2-1000. 

ENDIF 

Z2-ATAN ( -BSC2 ) 

ZC ( I )-0 . 5*(Z1+Z2 ) 
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AO  CONTINUE 

ZC(1)-ZC(2) 

ZC(M)-ZC(M-1) 

c* 

C*  DETERMINE  AVERAGE  CONTOUR  ORIENTATION  FOR  A  PROFILE  LINE 
C* 

KK-0 
ZCM-0 . 0 

DO  50  I-M-1,1,-1 

IF(D(I,JP).LT.O)  GOTO  99 
ZCM-ZCM+ZC ( I ) 

KK-KK+1 

50  CONTINUE 

99  ZCM— ZCM/REAL ( KK ) 

IF(ABS(ZCM) .LT.1.0E-6)  ZCM-0. 0 

RETURN 

END 
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c* 

C  ******'************************************************* 
SUBROUTINE  WAVEHD(KKK, D , J , H , E , V , DISS , ISTOP , NBR , 

&ZC.ZCM) 

C******************************************************* 

c* 

C*  THIS  ROUTINE  DETERMINES  THE  WAVE  HEIGHT  DISTRIBUTION 
C*  ACROSS-SHORE  FOR  A  SPECIFIC  PROFILE  J. 

C* 

INTEGER  MM,NN,NNGR,NM,M1 ,N1 ,NM1 , J , I 

PARAMETER  (MM-100 , NN-50 , NNGR-10 , NM-NN*MM , Ml-MM+1 , 

&N1-NN+1 , NM1— M1*N1) 

INTEGER  NBR (Ml ) ,M, ISTOP 

REAL  D(MM,NN) ,DHTOT(Ml) ,H(M1) , ZW(M1) , SXX(Ml) ,V(M1) , 

&SXY (Ml) , DMEAN(MM) , DH(M1) , DSEND , CC , PI , ZC(MM) , ZCC(MM) , 

&HBEG , ZBEG , T , LI , CGI , CN1 , E (Ml) , F(M1 ) , L2 , CG2 , CN2 , ZTEMP1 , 
&ZTEMP2 , GAMMA , KAPPA , AC , FSTAB , DI SS (MM) , DFS , ZCM , DCHG , LO 
LOGICAL  BREAK 

COMMON/BLOCKA/DX , DY , DTX , M , N , G , PI , GAMMA 
COMMON/BLOCKB/D50 ,  SRATIO ,  NY’ ,  SLOPE ,  DFS 
COMMON/BLOCKE/HBEG , ZBEG ,T , DCHG ,CC , DOFF 
DATA  KAPPA/0.15/ 

BREAK- . FALSE . 

DS END-0. 4*DFS 
C* 

C*  ZERO  ARRAY  CONTAINING  INFORMATION  ABOUT  BROKEN  WAVES 
C* 

DO  3  I-l.M+l 
NBR(I)— 0 
E(I)-0 

3  CONTINUE 
C* 

C*  DETERMINE  AVERAGE  DEPTH  BETWEEN  NEIGHBORING  PROFILES 
C*  AND  THE  LOCAL  CONTOUR  ORIENTATION 
C* 

CALL  CONPLA(D , J , DMEAN , ZCC , ZCM) 

DO  2  1-2, M 

DH( I )  — (DMEAN( I  - 1)+DMEAN(I ) )/2  . 

2  CONTINUE 

DH(M+1)— DMEAN(M) 

DH(l)-DMEANU) 

C* 

C*  CONVERT  LOCAL  COUNTOR  ANGLE  TO  REFERENCE  SYSTEM  WITH  RESPECT 
C*  TO  THE  X-AXIS  POINTING  OFFSHORE 
C* 

DO  5  I-l.M 

ZC(I)-PI/2 . -ZCC(I) 

5  CONTINUE 

ZCM-PI/2 . -ZCM 
C* 

C*  DETERMINE  WAVE  PROPERTIES  IN  THE  MOST  SEAWARD  CALCULATION  CELL 
C*  WHEN  THE  DEEPWATER  WAVE  CONDITIONS  ARE  GIVEN 
C* 

I F ( DOFF . LT . 0 ) THEN 
LO-1 . 5613*T**2 
FOFF— CC*HBEG**2*LO /T/2 
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ELSE 

CALL  LINWAV (Tt DOFF, LO , CGO , CNO) 

FOFF— CC*HBEG**2*CGO 
ENDIF 

CALL  LINWAV (T , DH(M+1) , LI , CGI , CN1) 

ZTEMPO-ZBEG-ZCM 
ZTEMP1-ASIN(L1/L0*SIN (ZTEMPO) ) 

ZW(M+1)-ZTEMP1+ZCM 

F (M+l) -FOFF*COS ( ZTEMPO ) /COS ( ZTEMP1 ) 

H(M+1)-SQRT(F(M+1)/CC/CG1) 

E(M+l)--PI*H(M+l)**2/4./Ll/SINH(4*PI*DH(M+l)/Ll) 
SXX(M+1)-CC*H(M+1)**2*(CN1*( (COS(ZW(M+l) -ZCM) )**2+l) -0 . 5) 
SXY(M+1 j-0 . 5*CC*HBEG**2*CN1*SIN( 2* (ZWCM+l ) - ZCM) ) 

C* 

C*  PROCEED  FROM  THE  END  OF  THE  GRID  IN  THE  SHOREWARD  DIRECTION 
C*  WHEN  CALCULATING  WAVE  PROPERTIES  ACROSS-SHORE 
C* 

DO  10  I-M.l.-l 
C* 

C*  DETERMINE  WAVE  PROPERTIES  FROM  LINEAR  WAVE  THEORY 
C* 

ETEST-E( 1+1 ) 

IF(DH(I) .LT.O)  ETEST-E(I+l)+E(I+l)-E(I+2) 

IF(DH( I )+ETEST . LT . 0)THEN 
ISTOP-I+1 
GOTO  77 
ENDIF 

CALL  LINWAV ( T , DH ( 1 ) +ETEST , L2 , CG2 , CN2 ) 

C* 

C*  CALCULATE  WAVE  REFRACTION  FROM  SNELLS  LAW  USING  THE  LOCAT. 

C*  CONTOUR  ORIENTATION 
C* 

ZTEMP1«ZW(I+1) -ZCM 
ZTEMP2-AS IN ( L2/L1*S IN ( ZTEMP1 ) ) 

ZW ( I )-ZTEMP2+ZCM 
C* 

C*  APPLY  RALLY'S  BREAKER  DECAY  MODEL  USING  THE  ENERGY  FLUX 
C* 

88  IF(H(I+1) .GT. GAMMA* (DH(I+1)+E(I+1) ) .OR. BREAK)THEN 

EMEAN-0 . 5*(ETEST+E(I+1) ) 

AC-KAPPA*DX/(DMEAN ( I ) +EMEAN ) 

FST AB-CC* ( 0 . 4* ( DMEAN ( I ) +EMEAN ) ) **  2* ( CG1+CG  2 ) / 2 . 
BREAK=.TRUE. 

NBR(I+1)-I 

C  WRITE(* , *)  'WAVE  ',1,'  IS  BROKEN' 

ELSE 

AC-0 

FSTAB-0 , 0 
NBR( 1+1 )— 0 
DISS(I)— 0 

C  WRITE(*,*)  'WAVE  ',1,'  IS  NOT  BROKEN’ 

ENDIF 


c* 

C*  CALCULATE  THE  ENERGY  FLUX 
C* 

F ( I)— (F(I+1)*(C0S (ZTEMP1 ) -0 . 5*AC)+AC*FSTAB)/(COS (ZTEMP2 ) 
&  +0 . 5*AC) 

IF(F(I) .LT.0)THEN 
ISTOP-I+1 
GOTO  77 
END  IF 
C* 

C*  CALCULATE  THE  ENERGY  DISSIPATION  PER  UNIT  VOLUME  IN  REGIONS  OF 
C*  BROKEN  WAVES 
C* 

IF (BREAK)  DISS(I)-KAPPA/(DMEAN(I)+E(I+1))**2*((F(I)+ 

&  F(I+l))/2. -FSTAB) 

C* 

C*  CHECK  IF  THE  WAVE  IS  REFORMING  BEFORE  REACHING  GRID  POINT  I 
C* 

IF((F(I)+F(I+l))/2. .LT. FSTAB) THEN 
BREAK-. FALSE. 

DISS(I)-0 

WRITE (*,*)  'SHOALING  -  HERE  WE  GO  AGAIN' 

GOTO  88 
END  IF 
C* 

C*  CALCULATE  THE  WAVE  HEIGHT  FROM  THE  ENERGY  FLUX 
C* 

H ( I ) -SQRT ( F ( I ) /CC/CG2 ) 

C* 

C*  CALCULATE  THE  RADIATION  STRESSES  SXX  AND  SXY 

C* 

SXX( I )-CC*H( I )**2*(CN2*( (COS (ZTEMP2) )**2+l ) -0 . 5) 
SXY(I)~0.5*CC*H(I)**2*CN2*SIN(2*ZTEMP2) 

C* 

C*  CALCULATE  THE  WAVE  SET-DOWN  OR  SET-UP 
C* 

E(I)-E(I+l)+(SXX(I+l)-SXX(I))/8./CC/(DMEAN(I)+E(I+l)) 
IF(DH(I)+E(I) . LT . 0 ) THEN 
ISTOP-I+1 
GOTO  77 
END  IF 
C* 

C*  UPDATE  WAVE  PARAMETERS 
C* 

CN1-CN2 

CG1-CG2 

L1-L2 

C  IF(I . GT . 1 . AND .DH(I-l).LT.O)  GOTO  99 

IF(DH(I)+E(I) .LT.DSEND)  GOTO  99 
10  CONTINUE 

C* 

C*  CALCULATE  THE  TOTAL  DEPTH 
C* 

99  IF(H(I) .GT. GAMMA*(DH( I )+E(l ) ) )  NBR(I)-1 

IF(NBR(I+1) . EQ.l. AND. 0 . 5*(F( I )+F( 1+1 ) ) -FSTAB. GT.O)  NBR(I)-1 
IF( I . GT . 0)THEN 
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ISTOP-I 

ELSE 

ISTOP-1 
END  IF 

77  DO  20  I-ISTOP.M+l 

DHTOT(I)-DH(I)+E(I) 

20  CONTINUE 

C* 

C*  WRITE  RESULT  TO  FILE 
C* 

IF(J .LE. 5 . AND .KKK. EQ . 1)THEN 

WRITE(24, 300)  'NR' , 'D' , ' DMEAN' , 'DM' , 'ZC' . 'ZCC' , 'ZW 
300  FORMAT (  '  '  , TR2  , A , TR4 ,  A,  TR4 ,  A , TR5 , A ,TR6 ,  A, TR5  ,  A , TR6  ,  A) 

WRITE(24 ,*) 

WRITE(21 ,400)  ’NR' , 'H’ , 'E' , 'SXX' , 'SXY' , 'DISS' , 'ZCM' 

400  FORMAT ( '  '  , TR2  ,  A ,TR4 , A , TR5 ,  A ,TR5  , A, TR4 ,  A,TR4 ,  A.TR4  ,  A) 
WRIT£(21 ,*) 

C  WRITE(22 ,*)  M-ISTOP+1 

CRD=180 . /PI 
DO  30  I=ISTOP,M 

Zl-90 . - ZC ( I ) *CRD 
Z2-90. -ZCC(I)*CRD 
Z3=90. -ZW(I)*CRD 

WRITE(24 , 100)  I , D ( I , J ) , DMEAN ( I ) , DH( I ) , Z1 , Z2 , Z3 
WRITE(21 , 200)  I,H(I) ,E(I) ,SXX(I) ,SXY(I) , DISS (I) , ZCM 
100  FORMAT ( '  ' ,I3,3F8.3,3F8.1) 

200  FORMAT ( '  ' ,I3,2F8.3,3F8.1,F8.3) 

C  WRITE(22 ,*)  REAL(I) ,H(I) 

30  CONTINUE 

ENDIF 

C  WRITE(22 , *)  M-ISTOP+1 

DO  40  I-ISTOP.M 

C  WRITE(22 ,*)  REAL(I) , -DH(I) 

40  CONTINUE 

C* 

C*  DETERMINE  THE  LONGSHORE  CURRENT 
CO¬ 
CALL  VELOC ( I STOP , NBR , H , ZW , ZC , ZCM , SXY , DHTOT , V ) 

RETURN 

END 
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c* 

c*******************************************,^.^.^.^^.^^ 
SUBROUTINE  LINWAV(T , DD , L,CG , CN) 
C******************************************************** 
c* 

C*  THIS  ROUTINE  CALCULATES  THE  WAVE  LENGTH,  GROUP  VELOCITY, 
C*  AND  THE  CONSTNANT  N  AT  A  SPECIFIED  DEPTH  D 
C* 

INTEGER  M,N 

REAL  T , DD , L , CG , CN , PI , G , DX , DY , DTX 
COMMON/BLOCKA/DX , DY , DTX , M , N , G , PI , GAMMA 
C* 

C*  USE  A  PADE  APPROXIMATION  TO  CALCULATE  THE  WAVE  LENGTH 
C* 

Y-(2*PI/T)**2*DD/G 

F“Y+1./(1 . +Y*(0. 66667+Y* (0. 3555+Y*(0. 16084+Y*(0. 06320+ 
&Y*(0 . 02174+Y*(0 . 00654+Y*(0 . 00171+Y*(0 . 00039+Y*0  000111 
&))))))))) 

L-2*PI/SQRT(Y*F/DD**2) 

C-4*PI*DD/L 

CN-0.5*(1+C/SINH(C)) 

CG=CN*L/T 

RETURN 

END 
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c* 

c ********************************************************* 
SUBROUTINE  VELOC(IS ,NBR,H, ZW, ZC , ZCM, SXY,DHTOT, V) 
C********************************************^^^^^^^ 

C* 

C*  CALCULATES  LONGSHORE  CURRENT  DISTRIBUTION  OVER  AN  ARBITRARY 
C*  BOTTOM  PROFILE;  WATER  DEPTH  &  WAVE  INFORMATION  GIVEN  EXTERNALLY 
C* 

INTEGER  I , I S , MM , NN , NNGR , NM , Ml , N1 , NM1 , M , N 
PARAMETER (MM-100 , NN-50 , NNGR-10 , NM-NN*MM 
&M1-MM+1 , Nl-NN+1 ,NM1-M1*N1) 

INTEGER  NBR(Ml) 

REAL  H(M1) , ZW(M1) , DHTOT (Ml) ,V(M1) ,UM(MI) , EPS (Ml) , EE(MM) 

&FF (MM) , DX , DTX , G , PI , AI , AHAT , BHAT , CHAT , RHAT , DEN , ZCM , DCHG 
&CAMMIX, CF, ALPHA, L,CG, CN, ZC(MM) ,HBEG,ZBEG  T  SXY(Ml)  KY 
COMMON/BLOCKA/DX , DY , DTX , M , N , G , PI , GAMMA 
COMMON/BLOCKE/HBEG , ZBEG , T , DCHG , CC , DOFF 
COMMON/BLOCKF/CF , GAMMIX , KY , QCUT 
ALPHA-0 . 0 
DXSQ=DX**2 
C* 

C*  BOTTOM  ORBITAL  VELOCITY 
C* 

DO  10  I  =  IS.M+l 

CALL  LINWAV (T , DHTOT(I ) , L, CG ,CN) 

UM( I )-0 . 5*H(I )*G*T/L/COSH ( 2*PI*DHT0T ( I ) /L) 

EPS ( I ) -GAMMIX*UM ( I ) *H ( I ) 

SXY ( I ) -SXY ( I ) /1000 . 

10  CONTINUE 
C* 

C*  SOLVE  FOR  VELOCITY 
C*  . 

C*  BOUNDARY  CONDITION’  FOR  VELOCITY  ON  SHOREWARD  SIDE  <  V(IS+1)  > 

C*  NOTE  X_AXIS  SHOULD  POINT  ONSHORE  FOR  DFN  OF  STRESSES 
C* 

ZI  =  0 . 5*(  ZW ( IS+1 )  +  ZW(IS))  -ZCM 
UMI-0 . 5* (UM( IS+1 )+UM(lS) ) 

V(IS)  =  - ( SXY (IS+1)  -  SXY ( IS) ) /DX 
&  /(2./PI*UMI*CF*(l  +  SIN(ZI  )**2) ) 

EE (IS)  -  0.0 
FF(IS)  -  V(IS) 

C* 

C*  FORWARD  SWEEP 
C* 

AI -EPS ( I S+ 1 ) *DHTOT ( I S+l ) 

DO  11  I  -  IS+1.  M 
IN-I+1 
IP-I-1 

UMI-0. 5*(UM(IN)+UM(I) ) 

AIN-EPS ( IN)*DHTOT (IN) 

AHAT  -  - AI/DXSQ 

CHAT  =  - AIN/DXSQ 

ZI  -  0 . 5*(ZW(IN)  +  ZW(I) ) -ZCM 

BI  -  2,/PI*CF*UMI*(l.  +  SIN(ZI)**2) 

BHAT  -  -(BI  +  (AI+AIN) /DXSQ) 

RHAT  -  (SXY (IN)  -  SXY ( I ) ) /DX  +  ALPHA 
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DEN  -  BHAT  -  AHAT*EE( IP) 

EE(I)  -  CHAT/DEN 

FF( I )  -  (AHAT*FF (IP)  +  RHAT)/DEN 
AI— AIN 
11  CONTINUE 
C* 

C*  BOUNDARY  CONDITION  FOR  VELOCITY  AT  SEAWARD  SIDE  <  V(M1)  > 
C* 

V(M+1)  -  0.0 

c* 

C*  BACKWARD  SWEEP 


C* 

DO  21  I  -  M,  IS+1,  -1 

V(I)  -  EE(I)*V(I+1)  +  FF(I) 

21  CONTINUE 
IVS-IVS+1 
C* 

C*  . . - .  END  OF  V  CALCULATION 

C* 

C  DO  55  I-IS.M 

C  WRITE(23, 100)  I , EPS (I) ,H(I) ,H(I+1) , UM( I ) , V( I ) 

C100  FORMAT ( '  ' , 13 , TR2 , F8 . 5 , 5 (TR2 , F8 . 3 ) ) 

C55  CONTINUE 

C  WRITE(23 ,*) 

C  WRITE(23 .*) 


RETURN 

END 
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c* 

£******'**■**•*■*■*■*'**■**•***■*  ******************************************* 
SUBROUTINE  TRANY ( I P , I STOP , NRU , NFS , NBR , V , H , ZC , ZCM , QY , 
&QYX) 

C*^*****************-*'*****'***'*'*'*-*********'*’******************-* 

c* 

C*  THIS  ROUTINE  DETERMINES  THE  LONGSHORE  TRANSPORT  RATE  ACROSS 
C*  THE  SURF  ZONE  ACCORDING  TO  THE  KRAUS  FORMULA. 

C* 

PARAMETER ( MM-lOO , NN-50 , NNGR-10 , NM-NN*MM , 

&M1-MM+1 .Nl-NN+l ,NM1-M1*N1) 

INTEGER  ISTOP,M,N,NBR(Ml) , IBRIGH, IBLEFT 
REAL  V(M1) ,H(M1) ,QY(M1,N1) , DX , DY , DTX.G , PI , KY , ZC(MM) , 
&QYX(M1,N1) 

COMMON  /BLOCKA/DX , DY , DTX , M , N , G , P I , GAMMA 
COMMON  /BLOCKF/CF , GAMMIX , KY , QCUT 
COMMON  /BLOCKG/IBRIGH, IBLEFT 
C* 

C*  DETERMINE  LONGSHORE  TRANSPORT  RATE 
C* 

IF(IP . GT , 1 . AND . IP . LT . N+l )THEN 
DO  5  I-l.ISTOP-l 
QY(I,IP)-0 
5  CONTINUE 

I F( I STOP. GT. NFS) THEN 

WRITE(* , *)  'TOO  BAD  -  NFS  <  ISTOP' 

ENDIF 

DO  10  I-ISTOP.M 
C  DO  10  I— NFS+5 ,M 

HAV-0.5*(H(I)+H(I+1)) 

QY ( I , I P ) "KY* V ( I ) *HAV - QCUT 
10  CONTINUE 

C  DO  33  I-NRU.NFS+4 

C  QY ( I , I P) -QY (NFS+5, IP) 

C33  CONTINUE 

ELSEIF(IP.EQ.1)THEN 

C* 

C*  SET  BOUNDARY  CONDITION  AT  RIGHT  END  OF  GRID 
C* 

DO  15  I-1,M 

IF( IBRIGH . EQ . 1 )  QY(I,l)-0 
15  CONTINUE 

ELSE 
C* 

C*  SET  BOUNDARY  CONDITION  AT  LEFT  END  OF  GRID 
C* 

DO  25  I-l.M 

I F( IBLEFT. EQ . 1)THEN 
QY(I,N+l)-0 
ELSE 

QY(I , N+l )“QY ( I , N) 

QYX ( I ,  N+ 1 ) -QYX ( I  ,N) 

ENDIF 

25  CONTINUE 

ENDIF 
C* 
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c*  MODIFY  TRANSPORT  RATE  ACCORDING  TO  LOCAL  CONTOUR 

IF(IBLEFT.NE.O.OR.IP.NE.N+l)THEN 

DO  35  I-l.M 

QY(I,IP)-<5Y(I1IP)*C0S(PI/2.  -ZCM) 

35  C0N?Se'IP)'QY<I,IP)*SIN(PI/2--ZCM) 

END  IF 

QYX(M+1 , IP)-0. 0 

RETURN 

END 
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APPENDIX  B 
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HEADING 

HEADING 

HEADING 

HEADING 

RUN  TITLE: 

2D-TEST  RUN 

B.  NUMBER  OF  CALCULATION  LOOPS  AND  CELLS:  NCL,  M(C-S),  N(L-S) 
160  100  10 

F.  TIME  AND  SPACE  INTERVAL:  DT(MIN) ,  DX(M).  DY(M) 

20  2  20 

LONGSHORE  TRANSPORT  PARAMETERS:  CF.GAMMIX 

0.01  1.0 

MEDIAN  SAND  GRAIN  SIZE:  D50(MM) 

0.2 

OUTPUT  INTERVAL  EXPRESSED  AS  NO.  OF  CALC.  TIME  LOOPS:  IOUT 
160 

N.  DEPTH  CORRESPONDING  TO  OFFSHORE  WAVE  DATA:  DZ 
2.96 

U.  SAND  TRANSPORT  CALIBRATION  COEFFICIENTS:  KY(L-S),  KX(C-S) 
5E-4  0.5 

BOUNDARY  CONDITIONS  (O-OPEN  COAST,  1-GROIN) 

1  0 
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